Module bianchi_sky_mod

module bianchi_sky_mod

        ! Uses
    use s2_types_mod
    use s2_sky_mod
    use s2_error_mod
    use bianchi_error_mod

        ! Types
    public type bianchi_sky

        ! Variables
    integer, public, parameter :: BIANCHI_SKY_QUAD_DIRECT = 1
    integer, public, parameter :: BIANCHI_SKY_QUAD_QTRAP = 2
    integer, public, parameter :: BIANCHI_SKY_QUAD_QSIMP = 3
    logical, private, parameter :: BIANCHI_DISABLE_PLM1TABLE = .false.
    real (kind=s2_dp), private, parameter :: BIANCHI_CMB_T = 1d0

        ! Subroutines and functions
    public function bianchi_sky_init (omega0, x, zE, s12H, s13H, rhand, nside, quad, N) result (b)
    public function bianchi_sky_init_alm (omega0, x, zE, s12H, s13H, rhand, lmax, quad_AB, quad_IAB, N, alpha, beta, gamma) result (b)
    public subroutine bianchi_sky_free (b)
    public subroutine bianchi_sky_param_write (b)
    public subroutine bianchi_sky_compute_map (b, nside)
    public subroutine bianchi_sky_compute_alm (b, lmax, mmax)
    public subroutine bianchi_sky_get_alm (b, alm)
    public subroutine bianchi_sky_apply_beam (b, fwhm, lmax)
    public subroutine bianchi_sky_write (b, filename, file_type, comment)
    public subroutine bianchi_sky_rotate (b, alpha, beta, gamma)
    private function bianchi_sky_comp_s (tau, tau0, theta0, h) result (s)
    private function bianchi_sky_comp_psi (tau, tau0, theta0, h) result (psi)
    private function bianchi_sky_comp_c1 (omega0, x) result (c1)
    private function bianchi_sky_comp_c2 (sE, zE) result (c2)
    private function bianchi_sky_comp_c3 (omega0, h) result (c3)
    private function bianchi_sky_comp_A (theta, tauE, tau0, h, zE, c1, c3, quad, N) result (Atheta)
    private function bianchi_sky_comp_B (theta, tauE, tau0, h, zE, c1, c3, quad, N) result (Btheta)
    private function bianchi_sky_comp_IA (tauE, tau0, h, zE, c1, c3, l, quad_AB, quad_IAB, N, A_grid) result (IA)
    private function bianchi_sky_comp_IB (tauE, tau0, h, zE, c1, c3, l, quad_AB, quad_IAB, N, B_grid) result (IB)
    private function bianchi_sky_integrand_A (tau, tau0, theta0, h) result (integrand)
    private function bianchi_sky_integrand_B (tau, tau0, theta0, h) result (integrand)
    private function bianchi_sky_integrand_IA (theta, tauE, tau0, h, zE, c1, c3, l, quad_AB, N, A_grid) result (integrand)
    private function bianchi_sky_integrand_IB (theta, tauE, tau0, h, zE, c1, c3, l, quad_AB, N, B_grid) result (integrand)
    private subroutine trapzd_a (func, tau0, theta0, h, a, b, s, n)
    private function qtrap_a (func, tau0, theta0, h, a, b) result (integral)
    private function qsimp_a (func, tau0, theta0, h, a, b) result (integral)
    private subroutine trapzd_ia (func, tauE, tau0, h, zE, c1, c3, l, quad_AB, a, b, s, n, N_quad)
    private function qtrap_ia (func, tauE, tau0, h, zE, c1, c3, l, quad_AB, a, b, N_quad) result (integral)
    private function qsimp_ia (func, tauE, tau0, h, zE, c1, c3, l, quad_AB, a, b, N_quad) result (integral)
    private function plgndr (l, m, x)
    private function asinh (x) result (y)

end module bianchi_sky_mod

Provides functionality to simulate a Bianchi VII_h model of the CMB. Uses the s2_sky module to create healpix sky maps.

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

Version: 0.1 June 2005

Description of Types

bianchi_sky

public type bianchi_sky
    private
    logical :: init = .false.
    real (kind=s2_dp) :: omega0 = 0.0d0
    real (kind=s2_dp) :: x = 0.0d0
    real (kind=s2_dp) :: zE = 0.0d0
    real (kind=s2_dp) :: s12H = 0.0d0
    real (kind=s2_dp) :: s13H = 0.0d0
    logical :: rhand = .true.
    real (kind=s2_dp) :: wH = 0.0d0
    real (kind=s2_dp) :: h = 0.0d0
    real (kind=s2_dp) :: tau0 = 0.0d0
    real (kind=s2_dp) :: tauE = 0.0d0
    type (s2_sky) :: sky
end type bianchi_sky

init: Initialisation status.
omega0: Density (input parameter).
x: Related to characteristic wavelength over which basis vectors change orientation (input parameter).
zE: Red shift (input parameter).
s12H: Normalised shear 12 (normalised to Hubble constant) (input parameter).
s13H: Normalised shear 13 (normalised to Hubble constant) (input parameter).
rhand: Logical specifying handedness of map (true=right).
wH: Normalised vorticity (normalised to Hubble constant).
h: Related to characteristic wavelength over which basis vectors change orientation.
tau0: Conformal time of photon reception.
tauE: Conformal time of photon emission.
sky: Simulated map.

Description of Variables

BIANCHI_SKY_QUAD_DIRECT

integer, public, parameter :: BIANCHI_SKY_QUAD_DIRECT = 1

Quadrature type: Direct (rectangular).

BIANCHI_SKY_QUAD_QTRAP

integer, public, parameter :: BIANCHI_SKY_QUAD_QTRAP = 2

Quadrature type: Trapezium rule.

BIANCHI_SKY_QUAD_QSIMP

integer, public, parameter :: BIANCHI_SKY_QUAD_QSIMP = 3

Quadrature type: Simpson's rule.

BIANCHI_DISABLE_PLM1TABLE

logical, private, parameter :: BIANCHI_DISABLE_PLM1TABLE = .false.

BIANCHI_CMB_T

real (kind=s2_dp), private, parameter :: BIANCHI_CMB_T = 1d0

Description of Subroutines and Functions

bianchi_sky_init

public function bianchi_sky_init (omega0, x, zE, s12H, s13H, rhand, nside, quad, N) result (b)
    real (kind=s2_dp), intent(in) :: omega0
    real (kind=s2_dp), intent(in) :: x
    real (kind=s2_dp), intent(in) :: zE
    real (kind=s2_dp), intent(in) :: s12H
    real (kind=s2_dp), intent(in) :: s13H
    logical, intent(in) :: rhand
    integer, intent(in) :: nside
    integer, intent(in) :: quad
    integer, optional, intent(in) :: N
    type (bianchi_sky) :: b
    ! Calls: bianchi_error, in_ring, pix2ang_ring
end function bianchi_sky_init

Initialise bianchi object by performing a bianchi simulation, computed in real space, with the specified parameters.

Variables:

omega0: Input omega0 parameter (see bianchi data type for explanation).
x: Input x parameter (see bianchi data type for explanation).
zE: Input zE parameter (see bianchi data type for explanation).
s12H: Input s12H parameter (see bianchi data type for explanation).
s13H: Input s13H parameter (see bianchi data type for explanation).
rhand: Logical to specify handedness of map.
nside: Nside of Healpix map to generate.
quad: Quadrature rule to use.
[N]: Resolution of numerical integration (number of terms used to evaluate integral). If not specified default of 10 is used.
b: Initialised bianchi object with all parameters and simulated map calculated.

Author: J. D. McEwen

Version: 0.1 June 2005

bianchi_sky_init_alm

public function bianchi_sky_init_alm (omega0, x, zE, s12H, s13H, rhand, lmax, quad_AB, quad_IAB, N, alpha, beta, gamma) result (b)
    real (kind=s2_dp), intent(in) :: omega0
    real (kind=s2_dp), intent(in) :: x
    real (kind=s2_dp), intent(in) :: zE
    real (kind=s2_dp), intent(in) :: s12H
    real (kind=s2_dp), intent(in) :: s13H
    logical, intent(in) :: rhand
    integer, intent(in) :: lmax
    integer, intent(in) :: quad_AB
    integer, intent(in) :: quad_IAB
    integer, optional, intent(in) :: N
    real (kind=s2_sp), optional, intent(in) :: alpha
    real (kind=s2_sp), optional, intent(in) :: beta
    real (kind=s2_sp), optional, intent(in) :: gamma
    type (bianchi_sky) :: b
    ! Calls: bianchi_error, s2_dl_beta_operator
end function bianchi_sky_init_alm

Initialise bianchi object by performing a bianchi simulation, computed in harmonic space, with the specified parameters.

Notes:

Uses lookup table for Plms for N=100, quad_IA=direct.

Variables:

omega0: Input omega0 parameter (see bianchi data type for explanation).
x: Input x parameter (see bianchi data type for explanation).
zE: Input zE parameter (see bianchi data type for explanation).
s12H: Input s12H parameter (see bianchi data type for explanation).
s13H: Input s13H parameter (see bianchi data type for explanation).
rhand: Logical to specify handedness of map.
nside: Nside of Healpix map to generate.
lmax: Maximum harmonic l to compute alms up to.
quad_AB: Quadrature rule to use for AB integration.
quad_IAB: Quadrature rule to use for IAB integration.
[N]: Resolution of numerical integration (number of terms used to evaluate integral). If not specified default of 10 is used.
[beam_fwhm]: Full-width-half-maximum of Gaussian beam to apply. If not present then no beam is applied. (Must be specified in arcmin.)
b: Initialised bianchi object with all parameters and simulated map calculated.

Author: J. D. McEwen

Version: 0.1 July 2005

bianchi_sky_free

public subroutine bianchi_sky_free (b)
    type (bianchi_sky), intent(inout) :: b
    ! Calls: bianchi_error, s2_sky_free
end subroutine bianchi_sky_free

Free all memory associated with a bianchi object.

Variables:

b: Bianchi object to free.

Author: J. D. McEwen

Version: 0.1 June 2005

bianchi_sky_param_write

public subroutine bianchi_sky_param_write (b)
    type (bianchi_sky), intent(in) :: b
end subroutine bianchi_sky_param_write

Write parameters of the Bianchi simulation to standard output.

Variables:

b: Bianchi object containing parameter atrtributes to write.

Author: J. D. McEwen

Version: 0.1 June 2005

bianchi_sky_compute_map

public subroutine bianchi_sky_compute_map (b, nside)
    type (bianchi_sky), intent(inout) :: b
    integer, intent(in) :: nside
    ! Calls: bianchi_error, s2_sky_compute_map
end subroutine bianchi_sky_compute_map

Compute the map of the bianchi sky, assuming the alms are already defined.

Variables:

b: Bianchi object containing alms to compute map of.
nside: Healpix nside to compute map at.

Author: J. D. McEwen

Version: 0.1 July 2005

bianchi_sky_compute_alm

public subroutine bianchi_sky_compute_alm (b, lmax, mmax)
    type (bianchi_sky), intent(inout) :: b
    integer, intent(in) :: lmax
    integer, intent(in) :: mmax
    ! Calls: bianchi_error, s2_sky_compute_alm
end subroutine bianchi_sky_compute_alm

Compute the alm of the bianchi sky, assuming the map is already defined.

Variables:

b: Bianchi object containing sky to compute alms of.
lmax: Maximum harmonic l to consider when computing alms.
mmax: Maximum harmonic m to consider when computing alms.

Author: J. D. McEwen

Version: 0.1 July 2005

bianchi_sky_get_alm

public subroutine bianchi_sky_get_alm (b, alm)
    type (bianchi_sky), intent(in) :: b
    complex (kind=s2_spc), intent(out), dimension (:,:) :: alm
    ! Calls: bianchi_error, s2_sky_get_alm
end subroutine bianchi_sky_get_alm

Get alms contained in a bianchi sky object. Alms must already be computed.

Notes:

Error occurs if alms are not already computed.

Variables:

b: Bianchi object containing sky that in turn contained the alms to get.

Author: J. D. McEwen

Version: 0.1 July 2005

bianchi_sky_apply_beam

public subroutine bianchi_sky_apply_beam (b, fwhm, lmax)
    type (bianchi_sky), intent(inout) :: b
    real (kind=s2_sp), intent(in) :: fwhm
    integer, intent(in) :: lmax
    ! Calls: s2_pl_free, s2_sky_conv
end subroutine bianchi_sky_apply_beam

Apply Gaussian beam with specified FWHM. (FWHM must be passed in arcmin.)

Notes:

Error occurs if alms are not already computed.

Variables:

b: Bianchi object containing sky that is to be comvolved with the beam.
fwhm: Gaussian beam FWHM to use (specified in arcmin).
lmax: Maximum harmonic l to consider.

Author: J. D. McEwen

Version: 0.1 July 2005

bianchi_sky_write

public subroutine bianchi_sky_write (b, filename, file_type, comment)
    type (bianchi_sky), intent(inout) :: b
    character (len=*), intent(in) :: filename
    integer, intent(in) :: file_type
    character (len=*), optional, intent(in) :: comment
    ! Calls: bianchi_error, s2_error, s2_sky_io_fits_write, s2_sky_write_map_file
end subroutine bianchi_sky_write

Write the bianchi simulated sky to a file.

Variables:

b: Bianchi object to save sky of.
filename: Name of output file.
file_type: Type of output file, either fits map or sky file (see s2_sky_mod for more details). Integer flag that may be either S2_SKY_FILE_TYPE_MAP or S2_SKY_FILE_TYPE_SKY.
[comment]: Optional comment to append to file header.

Author: J. D. McEwen

Version: 0.1 June 2005

bianchi_sky_rotate

public subroutine bianchi_sky_rotate (b, alpha, beta, gamma)
    type (bianchi_sky), intent(inout) :: b
    real (kind=s2_sp) :: alpha
    real (kind=s2_sp) :: beta
    real (kind=s2_sp) :: gamma
    ! Calls: bianchi_error, s2_sky_rotate
end subroutine bianchi_sky_rotate

Rotate the simulated sky of the bianchi object.

Variables:

b: Bianchi object containing sky to be rotated.
alpha: Alpha Euler angle of the rotation.
beta: Beta Euler angle of the rotation.
gamma: Gamma Euler angle of the rotation.

Author: J. D. McEwen

Version: 0.1 June 2005

bianchi_sky_comp_s

private function bianchi_sky_comp_s (tau, tau0, theta0, h) result (s)
    real (kind=s2_dp), intent(in) :: tau
    real (kind=s2_dp), intent(in) :: tau0
    real (kind=s2_dp), intent(in) :: theta0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp) :: s
end function bianchi_sky_comp_s

Compute s variable in bianchi sky simulation.

Variables:

tau: Conformal time to evaluate s for.
tau0: tau0 parameter (see bianchi data type for explanation).
theta0: theta0 position to evaluate s for.
h: h parameter (see bianchi data type for explanation).
s: s function value evaluated.

Author: J. D. McEwen

Version: 0.1 June 2005

bianchi_sky_comp_psi

private function bianchi_sky_comp_psi (tau, tau0, theta0, h) result (psi)
    real (kind=s2_dp), intent(in) :: tau
    real (kind=s2_dp), intent(in) :: tau0
    real (kind=s2_dp), intent(in) :: theta0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp) :: psi
end function bianchi_sky_comp_psi

Compute psi variable in bianchi sky simulation.

Variables:

tau: Conformal time to evaluate psi for.
tau0: tau0 parameter (see bianchi data type for explanation).
theta0: theta0 position to evaluate psi for.
h: h parameter (see bianchi data type for explanation).
psi: psi function value evaluated.

Author: J. D. McEwen

Version: 0.1 June 2005

bianchi_sky_comp_c1

private function bianchi_sky_comp_c1 (omega0, x) result (c1)
    real (kind=s2_dp), intent(in) :: omega0
    real (kind=s2_dp), intent(in) :: x
    real (kind=s2_dp) :: c1
end function bianchi_sky_comp_c1

Compute c1 variable in bianchi sky simulation.

Variables:

omega0: omega0 parameter (see bianchi data type for explanation).
x: x parameter (see bianchi data type for explanation).
c1: Value of constant c1 computed.

Author: J. D. McEwen

Version: 0.1 June 2005

bianchi_sky_comp_c2

private function bianchi_sky_comp_c2 (sE, zE) result (c2)
    real (kind=s2_dp), intent(in) :: sE
    real (kind=s2_dp), intent(in) :: zE
    real (kind=s2_dp) :: c2
end function bianchi_sky_comp_c2

Compute c2 variable in bianchi sky simulation.

Variables:

sE: sE parameter (value of s function computed at emission).
zE: zE parameter (see bianchi data type for explanation).
c2: Value of constant c2 computed.

Author: J. D. McEwen

Version: 0.1 June 2005

bianchi_sky_comp_c3

private function bianchi_sky_comp_c3 (omega0, h) result (c3)
    real (kind=s2_dp), intent(in) :: omega0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp) :: c3
end function bianchi_sky_comp_c3

Compute c3 variable in bianchi sky simulation.

Variables:

omega0: omega0 parameter (see bianchi data type for explanation).
h: h parameter (see bianchi data type for explanation).
c3: Value of constant c3 computed.

Author: J. D. McEwen

Version: 0.1 June 2005

bianchi_sky_comp_A

private function bianchi_sky_comp_A (theta, tauE, tau0, h, zE, c1, c3, quad, N) result (Atheta)
    real (kind=s2_dp), intent(in) :: theta
    real (kind=s2_dp), intent(in) :: tauE
    real (kind=s2_dp), intent(in) :: tau0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp), intent(in) :: zE
    real (kind=s2_dp), intent(in) :: c1
    real (kind=s2_dp), intent(in) :: c3
    integer, intent(in) :: quad
    integer, optional, intent(in) :: N
    real (kind=s2_dp) :: Atheta
    ! Calls: bianchi_error
end function bianchi_sky_comp_A

Compute A(theta) for Bianchi simulation.

Variables:

theta: theta value to eavluate B(theta) for.
- tauE: tauE parameter (see bianchi data type for explanation).
- tau0: tau0 parameter (see bianchi data type for explanation).
- h: h parameter (see bianchi data type for explanation).
- zE: zE parameter (see bianchi data type for explanation).
- c1: Bianchi simulation c1 parameter (constant parameter).
- c3: Bianchi simulation c3 parameter (constant parameter).
- quad: Integer specifying quadrature rule to use when evaluating integrals.
- [N]: Number of terms in use in direct quadrature.
- Atheta: Value of Atheta term evaluated.

Author: J. D. McEwen

Version: 0.1 July 2005

bianchi_sky_comp_B

private function bianchi_sky_comp_B (theta, tauE, tau0, h, zE, c1, c3, quad, N) result (Btheta)
    real (kind=s2_dp), intent(in) :: theta
    real (kind=s2_dp), intent(in) :: tauE
    real (kind=s2_dp), intent(in) :: tau0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp), intent(in) :: zE
    real (kind=s2_dp), intent(in) :: c1
    real (kind=s2_dp), intent(in) :: c3
    integer, intent(in) :: quad
    integer, optional, intent(in) :: N
    real (kind=s2_dp) :: Btheta
    ! Calls: bianchi_error
end function bianchi_sky_comp_B

Compute B(theta) for Bianchi simulation.

Variables:

theta: theta value to eavluate B(theta) for.
- tauE: tauE parameter (see bianchi data type for explanation).
- tau0: tau0 parameter (see bianchi data type for explanation).
- h: h parameter (see bianchi data type for explanation).
- zE: zE parameter (see bianchi data type for explanation).
- c1: Bianchi simulation c1 parameter (constant parameter).
- c3: Bianchi simulation c3 parameter (constant parameter).
- quad: Integer specifying quadrature rule to use when evaluating integrals.
- [N]: Number of terms in use in direct quadrature.
- Btheta: Value of Btheta term evaluated.

Author: J. D. McEwen

Version: 0.1 July 2005

bianchi_sky_comp_IA

private function bianchi_sky_comp_IA (tauE, tau0, h, zE, c1, c3, l, quad_AB, quad_IAB, N, A_grid) result (IA)
    real (kind=s2_dp), intent(in) :: tauE
    real (kind=s2_dp), intent(in) :: tau0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp), intent(in) :: zE
    real (kind=s2_dp), intent(in) :: c1
    real (kind=s2_dp), intent(in) :: c3
    integer, intent(in) :: l
    integer, intent(in) :: quad_AB
    integer, intent(in) :: quad_IAB
    integer, optional, intent(in) :: N
    real (kind=s2_dp), optional, intent(in), dimension (0:) :: A_grid
    real (kind=s2_dp) :: IA
    ! Calls: bianchi_error
end function bianchi_sky_comp_IA

Compute IA_l integral required in Bianchi simulation.

Variables:

tauE: tauE parameter (see bianchi data type for explanation).
tau0: tau0 parameter (see bianchi data type for explanation).
h: h parameter (see bianchi data type for explanation).
zE: zE parameter (see bianchi data type for explanation).
c1: Bianchi simulation c1 parameter (constant parameter).
c3: Bianchi simulation c3 parameter (constant parameter).
l: Harmonic l to compute IA_l for.
quad: Integer specifying quadrature rule to use when evaluating integrals (required to give to routine used to compute A(theta)).
[N]: Number of terms in use in direct quadrature.
IA: Value of IA_l integral evaluated.

Author: J. D. McEwen

Version: 0.1 July 2005

bianchi_sky_comp_IB

private function bianchi_sky_comp_IB (tauE, tau0, h, zE, c1, c3, l, quad_AB, quad_IAB, N, B_grid) result (IB)
    real (kind=s2_dp), intent(in) :: tauE
    real (kind=s2_dp), intent(in) :: tau0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp), intent(in) :: zE
    real (kind=s2_dp), intent(in) :: c1
    real (kind=s2_dp), intent(in) :: c3
    integer, intent(in) :: l
    integer, intent(in) :: quad_AB
    integer, intent(in) :: quad_IAB
    integer, optional, intent(in) :: N
    real (kind=s2_dp), optional, intent(in), dimension (0:) :: B_grid
    real (kind=s2_dp) :: IB
    ! Calls: bianchi_error
end function bianchi_sky_comp_IB

Compute IB_l integral required in Bianchi simulation.

Variables:

tauE: tauE parameter (see bianchi data type for explanation).
tau0: tau0 parameter (see bianchi data type for explanation).
h: h parameter (see bianchi data type for explanation).
zE: zE parameter (see bianchi data type for explanation).
c1: Bianchi simulation c1 parameter (constant parameter).
c3: Bianchi simulation c3 parameter (constant parameter).
l: Harmonic l to compute IB_l for.
quad: Integer specifying quadrature rule to use when evaluating integrals (required to give to routine used to compute B(theta)).
[N]: Number of terms in use in direct quadrature.
IB: Value of IB_l integral evaluated.

Author: J. D. McEwen

Version: 0.1 July 2005

bianchi_sky_integrand_A

private function bianchi_sky_integrand_A (tau, tau0, theta0, h) result (integrand)
    real (kind=s2_dp), intent(in) :: tau
    real (kind=s2_dp), intent(in) :: tau0
    real (kind=s2_dp), intent(in) :: theta0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp) :: integrand
end function bianchi_sky_integrand_A

Compute integrand of A integral.

Variables:

tau: Conformal time to evaluate integrand for.
tau0: tau0 parameter (see bianchi data type for explanation).
theta0: theta0 position to evaluate integrand for.
h: h parameter (see bianchi data type for explanation).
integrand: Integrand value evaluated.

Author: J. D. McEwen

Version: 0.1 June 2005

bianchi_sky_integrand_B

private function bianchi_sky_integrand_B (tau, tau0, theta0, h) result (integrand)
    real (kind=s2_dp), intent(in) :: tau
    real (kind=s2_dp), intent(in) :: tau0
    real (kind=s2_dp), intent(in) :: theta0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp) :: integrand
end function bianchi_sky_integrand_B

Compute integrand of B integral.

Variables:

tau: Conformal time to evaluate integrand for .
tau0: tau0 parameter (see bianchi data type for explanation).
theta0: theta0 position to evaluate integrand for.
h: h parameter (see bianchi data type for explanation).
integrand: Integrand value evaluated.

Author: J. D. McEwen

Version: 0.1 June 2005

bianchi_sky_integrand_IA

private function bianchi_sky_integrand_IA (theta, tauE, tau0, h, zE, c1, c3, l, quad_AB, N, A_grid) result (integrand)
    real (kind=s2_dp), intent(in) :: theta
    real (kind=s2_dp), intent(in) :: tauE
    real (kind=s2_dp), intent(in) :: tau0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp), intent(in) :: zE
    real (kind=s2_dp), intent(in) :: c1
    real (kind=s2_dp), intent(in) :: c3
    integer, intent(in) :: l
    integer, intent(in) :: quad_AB
    integer, optional, intent(in) :: N
    real (kind=s2_dp), optional, intent(in), dimension (0:) :: A_grid
    real (kind=s2_dp) :: integrand
    ! Calls: bianchi_error
end function bianchi_sky_integrand_IA

Compute integrand of IA_l integral.

Notes:

Computation of this integrand required the computation of A(theta) and the associated integral requried to compute A(theta).

Variables:

theta: theta angle to evaluate integrand for.
tauE: tauE parameter (see bianchi data type for explanation).
tau0: tau0 parameter (see bianchi data type for explanation).
h: h parameter (see bianchi data type for explanation).
zE: zE parameter (see bianchi data type for explanation).
c1: Bianchi simulation c1 parameter (constant parameter).
c3: Bianchi simulation c3 parameter (constant parameter).
l: Harmonic l to compute IA_l for.
quad: Integer specifying quadrature rule to use when evaluating integrals (required to give to routine used to compute A(theta)).
[N]: Number of terms in use in direct quadrature.
integrand: Integrand value evaluated.

Author: J. D. McEwen

Version: 0.1 July 2005

bianchi_sky_integrand_IB

private function bianchi_sky_integrand_IB (theta, tauE, tau0, h, zE, c1, c3, l, quad_AB, N, B_grid) result (integrand)
    real (kind=s2_dp), intent(in) :: theta
    real (kind=s2_dp), intent(in) :: tauE
    real (kind=s2_dp), intent(in) :: tau0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp), intent(in) :: zE
    real (kind=s2_dp), intent(in) :: c1
    real (kind=s2_dp), intent(in) :: c3
    integer, intent(in) :: l
    integer, intent(in) :: quad_AB
    integer, optional, intent(in) :: N
    real (kind=s2_dp), optional, intent(in), dimension (0:) :: B_grid
    real (kind=s2_dp) :: integrand
    ! Calls: bianchi_error
end function bianchi_sky_integrand_IB

Compute integrand of IB_l integral.

Notes:

Computation of this integrand required the computation of B(theta) and the associated integral requried to compute B(theta).

Variables:

theta: theta angle to evaluate integrand for.
tauE: tauE parameter (see bianchi data type for explanation).
tau0: tau0 parameter (see bianchi data type for explanation).
h: h parameter (see bianchi data type for explanation).
zE: zE parameter (see bianchi data type for explanation).
c1: Bianchi simulation c1 parameter (constant parameter).
c3: Bianchi simulation c3 parameter (constant parameter).
l: Harmonic l to compute IB_l for.
quad: Integer specifying quadrature rule to use when evaluating integrals (required to give to routine used to compute B(theta)).
[N]: Number of terms in use in direct quadrature.
integrand: Integrand value evaluated.

Author: J. D. McEwen

Version: 0.1 July 2005

trapzd_a

private subroutine trapzd_a (func, tau0, theta0, h, a, b, s, n)
    interface func
        function func (tau, tau0, theta0, h) result (integrand)
            real (kind=s2_dp), intent(in) :: tau
            real (kind=s2_dp), intent(in) :: tau0
            real (kind=s2_dp), intent(in) :: theta0
            real (kind=s2_dp), intent(in) :: h
            real (kind=s2_dp) :: integrand
        end function func
    end interface func
    real (kind=s2_dp), intent(in) :: tau0
    real (kind=s2_dp), intent(in) :: theta0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp), intent(IN) :: a
    real (kind=s2_dp), intent(IN) :: b
    real (kind=s2_dp), intent(INOUT) :: s
    integer, intent(IN) :: n
end subroutine trapzd_a

Computes nth stage of refinement of extended trapezoidal rule. Adapted from numerical recipes for the application at hand. For use in computing integrals required to compute A(theta) and B(theta).

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 A or B, (i.e. also takes parameters required to compute these integrals).
For our application tau0 and b will be the same, but clear and more general to program as two separate variables.

Variables:

func: "Pointer" to integrand function.
tau0: tau0 parameter (see bianchi data type for explanation).
theta0: theta0 position to evaluate integral refinement for.
h: h parameter (see bianchi data type for explanation).
a: Lower limit to evalue definite integral for.
b: Upper limit to evalue definite integral for.
s: Value of integral computed to current order.
n: Refinement order.

Author: J. D. McEwen

Version: 0.1 June 2005

qtrap_a

private function qtrap_a (func, tau0, theta0, h, a, b) result (integral)
    interface func
        function func (tau, tau0, theta0, h) result (integrand)
            real (kind=s2_dp), intent(in) :: tau
            real (kind=s2_dp), intent(in) :: tau0
            real (kind=s2_dp), intent(in) :: theta0
            real (kind=s2_dp), intent(in) :: h
            real (kind=s2_dp) :: integrand
        end function func
    end interface func
    real (kind=s2_dp), intent(in) :: tau0
    real (kind=s2_dp), intent(in) :: theta0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp), intent(IN) :: a
    real (kind=s2_dp), intent(IN) :: b
    real (kind=s2_dp) :: integral
    ! Calls: bianchi_error, trapzd_a
end function qtrap_a

Computes the integral of the function func from a to b using the trapezoid method. Adapted from numerical recipes for the application at hand. For use in computing integrals required to compute A(theta) and B(theta).

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 A or B, (i.e. also takes parameters required to compute these integrals).
For our application tau0 and b will be the same, but clear and more general to program as two separate variables.!!

Variables:

func: "Pointer" to integrand function.
tau0: tau0 parameter (see bianchi data type for explanation).
theta0: theta0 position to evaluate integral refinement for.
h: h parameter (see bianchi data type for explanation).
a: Lower limit to evalue definite integral for.
b: Upper limit to evalue definite integral for.
integral: Value of the evaluated integral.

Author: J. D. McEwen

Version: 0.1 June 2005

qsimp_a

private function qsimp_a (func, tau0, theta0, h, a, b) result (integral)
    interface func
        function func (tau, tau0, theta0, h) result (integrand)
            real (kind=s2_dp), intent(in) :: tau
            real (kind=s2_dp), intent(in) :: tau0
            real (kind=s2_dp), intent(in) :: theta0
            real (kind=s2_dp), intent(in) :: h
            real (kind=s2_dp) :: integrand
        end function func
    end interface func
    real (kind=s2_dp), intent(in) :: tau0
    real (kind=s2_dp), intent(in) :: theta0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp), intent(IN) :: a
    real (kind=s2_dp), intent(IN) :: b
    real (kind=s2_dp) :: integral
    ! Calls: bianchi_error, trapzd_a
end function qsimp_a

Computes the integral of the function func from a to b using Simpson's rule. Adapted from numerical recipes for the application at hand. For use in computing integrals required to compute A(theta) and B(theta).

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 A or B, (i.e. also takes parameters required to compute these integrals).
For our application tau0 and b will be the same, but clear and more general to program as two separate variables.

Variables:

func: "Pointer" to integrand function.
tau0: tau0 parameter (see bianchi data type for explanation).
theta0: theta0 position to evaluate integral refinement for.
h: h parameter (see bianchi data type for explanation).
a: Lower limit to evalue definite integral for.
b: Upper limit to evalue definite integral for.
integral: Value of the evaluated integral.

Author: J. D. McEwen

Version: 0.1 June 2005

trapzd_ia

private subroutine trapzd_ia (func, tauE, tau0, h, zE, c1, c3, l, quad_AB, a, b, s, n, N_quad)
    interface func
        function func (theta, tauE, tau0, h, zE, c1, c3, l, quad, N, AB_grid) result (integrand)
            real (kind=s2_dp), intent(in) :: theta
            real (kind=s2_dp), intent(in) :: tauE
            real (kind=s2_dp), intent(in) :: tau0
            real (kind=s2_dp), intent(in) :: h
            real (kind=s2_dp), intent(in) :: zE
            real (kind=s2_dp), intent(in) :: c1
            real (kind=s2_dp), intent(in) :: c3
            integer, intent(in) :: l
            integer, intent(in) :: quad
            integer, optional, intent(in) :: N
            real (kind=s2_dp), optional, intent(in), dimension (0:) :: AB_grid
            real (kind=s2_dp) :: integrand
        end function func
    end interface func
    real (kind=s2_dp), intent(in) :: tauE
    real (kind=s2_dp), intent(in) :: tau0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp), intent(in) :: zE
    real (kind=s2_dp), intent(in) :: c1
    real (kind=s2_dp), intent(in) :: c3
    integer, intent(in) :: l
    integer, intent(in) :: quad_AB
    real (kind=s2_dp), intent(IN) :: a
    real (kind=s2_dp), intent(IN) :: b
    real (kind=s2_dp), intent(INOUT) :: s
    integer, intent(IN) :: n
    integer, optional, intent(in) :: N_quad
end subroutine trapzd_ia

Computes nth stage of refinement of extended trapezoidal rule. Adapted from numerical recipes for the application at hand. For use in computing integrals IA_l and IB_l (although function name is only give the a subscript).

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 IA_l or IB_l, (i.e. also takes parameters required to compute these integrals).

Variables:

func: "Pointer" to integrand function.
tauE: tauE parameter (see bianchi data type for explanation).
tau0: tau0 parameter (see bianchi data type for explanation).
h: h parameter (see bianchi data type for explanation).
zE: zE parameter (see bianchi data type for explanation).
c1: Bianchi simulation c1 parameter (constant parameter).
c3: Bianchi simulation c3 parameter (constant parameter).
l: Harmonic l to compute integral for.
quad: Integer specifying quadrature rule to use when evaluating integral.
a: Lower limit to evalue definite integral for.
b: Upper limit to evalue definite integral for.
s: Value of integral computed to current order.
n: Refinement order.
[N_quad]: Number of terms in use in direct quadrature.

Author: J. D. McEwen

Version: 0.1 July 2005

qtrap_ia

private function qtrap_ia (func, tauE, tau0, h, zE, c1, c3, l, quad_AB, a, b, N_quad) result (integral)
    interface func
        function func (theta, tauE, tau0, h, zE, c1, c3, l, quad, N, AB_grid) result (integrand)
            real (kind=s2_dp), intent(in) :: theta
            real (kind=s2_dp), intent(in) :: tauE
            real (kind=s2_dp), intent(in) :: tau0
            real (kind=s2_dp), intent(in) :: h
            real (kind=s2_dp), intent(in) :: zE
            real (kind=s2_dp), intent(in) :: c1
            real (kind=s2_dp), intent(in) :: c3
            integer, intent(in) :: l
            integer, intent(in) :: quad
            integer, optional, intent(in) :: N
            real (kind=s2_dp), optional, intent(in), dimension (0:) :: AB_grid
            real (kind=s2_dp) :: integrand
        end function func
    end interface func
    real (kind=s2_dp), intent(in) :: tauE
    real (kind=s2_dp), intent(in) :: tau0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp), intent(in) :: zE
    real (kind=s2_dp), intent(in) :: c1
    real (kind=s2_dp), intent(in) :: c3
    integer, intent(in) :: l
    integer, intent(in) :: quad_AB
    real (kind=s2_dp), intent(IN) :: a
    real (kind=s2_dp), intent(IN) :: b
    integer, optional, intent(in) :: N_quad
    real (kind=s2_dp) :: integral
    ! Calls: bianchi_error, trapzd_ia
end function qtrap_ia

Computes the integral of the function func from a to b using the trapezoid method. Adapted from numerical recipes for the application at hand. For use in computing integrals IA_l and IB_l (although function name is only give the a subscript).

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 IA_l or IB_l, (i.e. also takes parameters required to compute these integrals).

Variables:

func: "Pointer" to integrand function.
tauE: tauE parameter (see bianchi data type for explanation).
tau0: tau0 parameter (see bianchi data type for explanation).
h: h parameter (see bianchi data type for explanation).
zE: zE parameter (see bianchi data type for explanation).
c1: Bianchi simulation c1 parameter (constant parameter).
c3: Bianchi simulation c3 parameter (constant parameter).
l: Harmonic l to compute integral for.
quad: Integer specifying quadrature rule to use when evaluating integral.
[N_quad]: Number of terms in use in direct quadrature.
a: Lower limit to evalue definite integral for.
b: Upper limit to evalue definite integral for.
integral: Value of the evaluated integral.

Author: J. D. McEwen

Version: 0.1 July 2005

qsimp_ia

private function qsimp_ia (func, tauE, tau0, h, zE, c1, c3, l, quad_AB, a, b, N_quad) result (integral)
    interface func
        function func (theta, tauE, tau0, h, zE, c1, c3, l, quad, N, AB_grid) result (integrand)
            real (kind=s2_dp), intent(in) :: theta
            real (kind=s2_dp), intent(in) :: tauE
            real (kind=s2_dp), intent(in) :: tau0
            real (kind=s2_dp), intent(in) :: h
            real (kind=s2_dp), intent(in) :: zE
            real (kind=s2_dp), intent(in) :: c1
            real (kind=s2_dp), intent(in) :: c3
            integer, intent(in) :: l
            integer, intent(in) :: quad
            integer, optional, intent(in) :: N
            real (kind=s2_dp), optional, intent(in), dimension (0:) :: AB_grid
            real (kind=s2_dp) :: integrand
        end function func
    end interface func
    real (kind=s2_dp), intent(in) :: tauE
    real (kind=s2_dp), intent(in) :: tau0
    real (kind=s2_dp), intent(in) :: h
    real (kind=s2_dp), intent(in) :: zE
    real (kind=s2_dp), intent(in) :: c1
    real (kind=s2_dp), intent(in) :: c3
    integer, intent(in) :: l
    integer, intent(in) :: quad_AB
    real (kind=s2_dp), intent(IN) :: a
    real (kind=s2_dp), intent(IN) :: b
    integer, optional, intent(in) :: N_quad
    real (kind=s2_dp) :: integral
    ! Calls: bianchi_error, trapzd_ia
end function qsimp_ia

Computes the integral of the function func from a to b using Simpson's rule. Adapted from numerical recipes for the application at hand. For use in computing integrals IA_l and IB_l (although function name is only give the a subscript).

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 A or B, (i.e. also takes parameters required to compute these integrals).

Variables:

func: "Pointer" to integrand function.
tauE: tauE parameter (see bianchi data type for explanation).
tau0: tau0 parameter (see bianchi data type for explanation).
h: h parameter (see bianchi data type for explanation).
zE: zE parameter (see bianchi data type for explanation).
c1: Bianchi simulation c1 parameter (constant parameter).
c3: Bianchi simulation c3 parameter (constant parameter).
l: Harmonic l to compute integral for.
quad: Integer specifying quadrature rule to use when evaluating integral.
[N_quad]: Number of terms in use in direct quadrature.
a: Lower limit to evalue definite integral for.
b: Upper limit to evalue definite integral for.
integral: Value of the evaluated integral.

Author: J. D. McEwen

Version: 0.1 June 2005

plgndr

private function plgndr (l, m, x)
    integer :: l
    integer :: m
    real (kind=s2_dp) :: x
    real (kind=s2_dp) :: plgndr
end function plgndr

Computes the associated Legendre function for l and m. Adapted from numerical recipes

Notes:

Numerical recipies comment: Computes the associated Legendre polynomial P_m^l(x).

Variables:

l: Legendre function l parameter.
m: Legendre function m parameter.
x: Point to evaluate specified Legendre funtion at.

Author: J. D. McEwen

Version: 0.1 June 2005

asinh

private function asinh (x) result (y)
    real (kind=s2_dp), intent(in) :: x
    real (kind=s2_dp) :: y
end function asinh

Implementation of asinh, since no intrinsic Fortran function.

Variables:

x: Asinh argument.
y: Asinh result.

Author: J. D. McEwen

Version: 0.1 June 2005