ssht_dl.c File Reference

#include <stdlib.h>
#include <math.h>
#include "ssht/ssht_types.h"
#include "ssht/ssht_error.h"
#include "ssht/ssht_dl.h"

Include dependency graph for ssht_dl.c:

Functions
double	logfact (int n)
double *	ssht_dl_calloc (int L, ssht_dl_size_t dl_size)
int	ssht_dl_get_offset (int L, ssht_dl_size_t dl_size)
int	ssht_dl_get_stride (int L, ssht_dl_size_t dl_size)
void	ssht_dl_beta_risbo_full_table (double *dl, double beta, int L, ssht_dl_size_t dl_size, int el, double *sqrt_tbl)
void	ssht_dl_beta_risbo_half_table (double *dl, double beta, int L, ssht_dl_size_t dl_size, int el, double *sqrt_tbl, double *signs)
void	ssht_dl_beta_risbo_eighth_table (double *dl, double beta, int L, ssht_dl_size_t dl_size, int el, double *sqrt_tbl, double *signs)
void	ssht_dl_beta_risbo_fill_eighth2quarter_table (double *dl, double *dl8, int L, ssht_dl_size_t dl_size, ssht_dl_size_t dl8_size, int el, double *signs)
void	ssht_dl_beta_kostelec_full_table (double *dlm1p1, double *dl, double beta, int L, ssht_dl_size_t dl_size, int el, double *sqrt_tbl, double *signs)
void	ssht_dl_beta_kostelec_line_table (double *dlm1p1_line, double *dl_line, double beta, int L, int mm, int el, double *sqrt_tbl, double *signs)
void	ssht_dl_beta_kostelec_halfline_table (double *dlm1p1_line, double *dl_line, double beta, int L, int mm, int el, double *sqrt_tbl, double *signs)
void	ssht_dl_halfpi_trapani_eighth_table (double *dl, int L, ssht_dl_size_t dl_size, int el, double *sqrt_tbl)
void	ssht_dl_halfpi_trapani_quarter_table (double *dl, int L, ssht_dl_size_t dl_size, int el, double *sqrt_tbl)
void	ssht_dl_halfpi_trapani_fill_eighth2righthalf_table (double *dl, int L, ssht_dl_size_t dl_size, int el, double *signs)
void	ssht_dl_halfpi_trapani_fill_eighth2quarter_table (double *dl, int L, ssht_dl_size_t dl_size, int el, double *signs)

Detailed Description

Functionality to compute Wigner plane and manage its storage.

Note

To access the allocated memory correctly, attention must be paid to its size. This depends of the size of the dl plane required, as specified by the ssht_dl_size_t type. Memory is allocated and accessed as follows for the following sizes.

• SSHT_DL_QUARTER:
  • allocated by dl = (double*)calloc(L*L, sizeof(double))
  • accessed by dl[m*L + mm], i.e. m in [0,(L-1)] and mm in [0,(L-1)]
• SSHT_DL_QUARTER_EXTENDED:
  • allocated by dl = (double*)calloc((L+1)*(L+3), sizeof(double))
  • require extended size when using Risbo recursion for partial Wigner plane
  • accessed by dl[(m + (L-1))*(L+3) + mm + L-1], i.e. m in [-(L-1),0] and mm in [-(L-1),0]
• SSHT_DL_HALF:
  • allocated by dl = (double*)calloc(L*(2*L-1), sizeof(double))
  • accessed by dl[m*(2*L-1) + mm + L-1], i.e. m in [0,(L-1)] and mm in [-(L-1),(L-1)]
• SSHT_DL_FULL:
  • allocated by dl = (double*)calloc((2*L-1)*(2*L-1), sizeof(double))
  • accessed by dl[(m + (L-1))*(2*L-1) + mm + L-1], i.e. m in [-(L-1),(L-1)] and mm in [-(L-1),(L-1)]

The routine ssht_dl_calloc is provided to allocate space to store a dl plane. The routines ssht_dl_get_offset and ssht_dl_get_stride are provided to get appropriate offsets and strides.

Author

Jason McEwen

Function Documentation

logfact()

double logfact

(

int

n

)

Computes the natural logarithm of an (integer) factorial.

Parameters

[in]

n

Integer to compute factorial of.

Return values

logfactn

Natural logarithm of factorial value computed.

Author

Numerical recipes.

ssht_dl_beta_kostelec_full_table()

void ssht_dl_beta_kostelec_full_table	(	double *	dlm1p1,
		double *	dl,
		double	beta,
		int	L,
		ssht_dl_size_t	dl_size,
		int	el,
		double *	sqrt_tbl,
		double *	signs
	)

Calculates (for m = -l:l and mm = -l:l) lth plane of a d-matrix for argument beta using the recursion method given in Kostelec and Rockmore (2010) (see equations (4.5)-(4.9)). For l>1, require the dl plane to be computed already with values for l-1 and l-2. Also takes a table of precomputed square roots of integers and signs to avoid recomputing them.

Parameters

[in,out]	dlm1p1	Wigner plane. On input this should be initialised to the plane computed for el-2. On output this will be replaced with the computed plane for el.
[in]	dl	Wigner plane already computed for el-1.
[in]	beta	Angle to compute Wigner plan for.
[in]	L	Harmonic band-limit.
[in]	dl_size	Size type of the memory to allocate.
[in]	el	Harmonic index to compute Wigner plane for.
[in]	sqrt_tbl	Precomputed array of square roots. The table element at index i should contain the value sqrt(i). Values from 0 to 2*el must be precomputed (i.e. sqrt_tbl should contian 2*el+1 elements).
[in]	signs	Precomputed array of signs. The array element at index i should contain the value (-1)^i. Values from 0 to el must be precomputed (i.e. signs should contian el+1 elements).

Return values

none

Author

Jason McEwen

ssht_dl_beta_kostelec_halfline_table()

void ssht_dl_beta_kostelec_halfline_table	(	double *	dlm1p1_line,
		double *	dl_line,
		double	beta,
		int	L,
		int	mm,
		int	el,
		double *	sqrt_tbl,
		double *	signs
	)

Calculates half line of d-matrix for all m = 0:l and given mm = -l:l for argument beta using the recursion method given in Kostelec and Rockmore (2010) (see equations (4.5)-(4.9)). For l>abs(mm), require the dl plane to be computed already with values for l-1 and l-2. Also takes a table of precomputed square roots of integers and signs to avoid recomputing them.

Parameters

[in,out]	dlm1p1_line	Wigner line. On input this should be initialised to the line computed for el-2. On output this will be replaced with the computed line for el.
[in]	dl	Wigner plane already computed for el-1.
[in]	beta	Angle to compute Wigner line for.
[in]	L	Harmonic band-limit.
[in]	mm	Azimuthal harmonic index to compute Wigner line for.
[in]	el	Harmonic index to compute Wigner plane for.
[in]	sqrt_tbl	Precomputed array of square roots. The table element at index i should contain the value sqrt(i). Values from 0 to 2*el must be precomputed (i.e. sqrt_tbl should contian 2*el+1 elements).
[in]	signs	Precomputed array of signs. The array element at index i should contain the value (-1)^i. Values from 0 to el must be precomputed (i.e. signs should contian el+1 elements).

Return values

none

Author

Jason McEwen

ssht_dl_beta_kostelec_line_table()

void ssht_dl_beta_kostelec_line_table	(	double *	dlm1p1_line,
		double *	dl_line,
		double	beta,
		int	L,
		int	mm,
		int	el,
		double *	sqrt_tbl,
		double *	signs
	)

Calculates line of d-matrix for all m = -l:l and given mm = -l:l for argument beta using the recursion method given in Kostelec and Rockmore (2010) (see equations (4.5)-(4.9)). For l>abs(mm), require the dl plane to be computed already with values for l-1 and l-2. Also takes a table of precomputed square roots of integers and signs to avoid recomputing them.

Parameters

[in,out]	dlm1p1_line	Wigner line. On input this should be initialised to the line computed for el-2. On output this will be replaced with the computed line for el.
[in]	dl	Wigner plane already computed for el-1.
[in]	beta	Angle to compute Wigner line for.
[in]	L	Harmonic band-limit.
[in]	mm	Azimuthal harmonic index to compute Wigner line for.
[in]	el	Harmonic index to compute Wigner plane for.
[in]	sqrt_tbl	Precomputed array of square roots. The table element at index i should contain the value sqrt(i). Values from 0 to 2*el must be precomputed (i.e. sqrt_tbl should contian 2*el+1 elements).
[in]	signs	Precomputed array of signs. The array element at index i should contain the value (-1)^i. Values from 0 to el must be precomputed (i.e. signs should contian el+1 elements).

Return values

none

Author

Jason McEwen

ssht_dl_beta_risbo_eighth_table()

void ssht_dl_beta_risbo_eighth_table	(	double *	dl,
		double	beta,
		int	L,
		ssht_dl_size_t	dl_size,
		int	el,
		double *	sqrt_tbl,
		double *	signs
	)

Calculates (for m = -l:0 and mm = -l:m) lth plane of a d-matrix for argument beta using Risbo's recursion method. For l>0, require the dl plane to be computed already with values for l-1. Also takes a table of precomputed square roots of integers and signs to avoid recomputing them.

Note

Risbo's recursion requires the boundaries of the eighth of the Wigner plane to be extended slightly. Consequently, this routine should only be used with dl_size either SSHT_DL_QUARTER_EXTENDED or SSHT_DL_FULL.

Parameters

[in,out]	dl	Wigner plane. On input this should be initialised to the plane computed for el-1. On output this will be replaced with the computed plane for el.
[in]	beta	Angle to compute Wigner line for.
[in]	L	Harmonic band-limit.
[in]	dl_size	Size type of the memory to allocate.
[in]	el	Harmonic index to compute Wigner plane for.
[in]	sqrt_tbl	Precomputed array of square roots. The table element at index i should contain the value sqrt(i). Values from 0 to 2*el must be precomputed (i.e. sqrt_tbl should contian 2*el+1 elements).
[in]	signs	Precomputed array of signs. The array element at index i should contain the value (-1)^i. Values from 0 to el must be precomputed (i.e. signs should contian el+1 elements).

Return values

none

Author

Jason McEwen

ssht_dl_beta_risbo_fill_eighth2quarter_table()

void ssht_dl_beta_risbo_fill_eighth2quarter_table	(	double *	dl,
		double *	dl8,
		int	L,
		ssht_dl_size_t	dl_size,
		ssht_dl_size_t	dl8_size,
		int	el,
		double *	signs
	)

Fill in quarter Wigner plane for m = 0:l and mm = 0:l from the eighth m = -l:0 and mm = -l:m. Takes a table of precomputed signs to avoid recomputing them.

Note

It is necessary to fill the quarter Wigner plane in separate memory since the eighth plane for m = -l:0 and mm = -l:m with be required in subsequent recursions.

Parameters

[out]	dl	Quarter Wigner plane, computed by copying dl8 and extending with symmetries.
[in]	dl8	Eighth Wigner plane.
[in]	L	Harmonic band-limit.
[in]	dl_size	Size type of the dl memory.
[in]	dl8_size	Size type of the dl8 memory.
[in]	el	Harmonic index to compute Wigner plane for.
[in]	signs	Precomputed array of signs. The array element at index i should contain the value (-1)^i. Values from 0 to el must be precomputed (i.e. signs should contian el+1 elements).

Return values

none

Author

Jason McEwen

ssht_dl_beta_risbo_full_table()

void ssht_dl_beta_risbo_full_table	(	double *	dl,
		double	beta,
		int	L,
		ssht_dl_size_t	dl_size,
		int	el,
		double *	sqrt_tbl
	)

Calculates (for m = -l:l and mm = -l:l) lth plane of a d-matrix for argument beta using Risbo's recursion method. For l>0, require the dl plane to be computed already with values for l-1. Also takes a table of precomputed square roots of integers to avoid recomputing them.

Parameters

[in,out]	dl	Wigner plane. On input this should be initialised to the plane computed for el-1. On output this will be replaced with the computed plane for el.
[in]	beta	Angle to compute Wigner line for.
[in]	L	Harmonic band-limit.
[in]	dl_size	Size type of the memory to allocate.
[in]	el	Harmonic index to compute Wigner plane for.
[in]	sqrt_tbl	Precomputed array of square roots. The table element at index i should contain the value sqrt(i). Values from 0 to 2*el must be precomputed (i.e. sqrt_tbl should contian 2*el+1 elements).

Return values

none

Author

Jason McEwen

ssht_dl_beta_risbo_half_table()

void ssht_dl_beta_risbo_half_table	(	double *	dl,
		double	beta,
		int	L,
		ssht_dl_size_t	dl_size,
		int	el,
		double *	sqrt_tbl,
		double *	signs
	)

Calculates (for m = -l:l and mm = -l:l) lth plane of a d-matrix for argument beta using Risbo's recursion method. Only half of the plane is computed by recusion and symmetry is used to fill the remaining plane. For l>0, require the dl plane to be computed already with values for l-1. Also takes a table of precomputed square roots of integers and signs to avoid recomputing them.

Parameters

[in,out]	dl	Wigner plane. On input this should be initialised to the plane computed for el-1. On output this will be replaced with the computed plane for el.
[in]	beta	Angle to compute Wigner line for.
[in]	L	Harmonic band-limit.
[in]	dl_size	Size type of the memory to allocate.
[in]	el	Harmonic index to compute Wigner plane for.
[in]	sqrt_tbl	Precomputed array of square roots. The table element at index i should contain the value sqrt(i). Values from 0 to 2*el must be precomputed (i.e. sqrt_tbl should contian 2*el+1 elements).
[in]	signs	Precomputed array of signs. The array element at index i should contain the value (-1)^i. Values from 0 to el must be precomputed (i.e. signs should contian el+1 elements).

Return values

none

Author

Jason McEwen

ssht_dl_calloc()

double* ssht_dl_calloc	(	int	L,
		ssht_dl_size_t	dl_size
	)

Allocate memory to store dl plane.

Note

See file note regarding storage and access of dl memory.

Parameters

[in]	L	Harmonic band-limit.
[in]	dl_size	Size type of the memory to allocate (specifying whether will need half or quarted dl plane, for example).

Return values

dl	Pointer to dl array with memory allocated.

Author

Jason McEwen

ssht_dl_get_offset()

int ssht_dl_get_offset	(	int	L,
		ssht_dl_size_t	dl_size
	)

Compute the offset used for accessing dl memory.

Note

See file note regarding storage and access of dl memory.

Parameters

[in]	L	Harmonic band-limit.
[in]	dl_size	Size type of the memory to allocate (specifying whether will need half or quarted dl plane, for example).

Return values

offset

The offset to be used to access dl memory.

Author

Jason McEwen

ssht_dl_get_stride()

int ssht_dl_get_stride	(	int	L,
		ssht_dl_size_t	dl_size
	)

Compute the stride used for accessing dl memory.

Note

See file note regarding storage and access of dl memory.

Parameters

[in]	L	Harmonic band-limit.
[in]	dl_size	Size type of the memory to allocate (specifying whether will need half or quarted dl plane, for example).

Return values

stride

The stride to be used to access dl memory.

Author

Jason McEwen

ssht_dl_halfpi_trapani_eighth_table()

void ssht_dl_halfpi_trapani_eighth_table	(	double *	dl,
		int	L,
		ssht_dl_size_t	dl_size,
		int	el,
		double *	sqrt_tbl
	)

Calculates eighth (for m = 0:l and mm = 0:m) of lth plane of a d-matrix for PI/2 using Trapani & Navaza's recursion method. For l>0, require the dl plane to be computed already with values for l-1. Also takes a table of precomputed square roots of integers to avoid recomputing them.

Parameters

[in,out]	dl	Wigner plane. On input this should be initialised to the plane computed for el-1. On output this will be replaced with the computed plane for el.
[in]	L	Harmonic band-limit.
[in]	dl_size	Size type of the memory to allocate.
[in]	el	Harmonic index to compute Wigner plane for.
[in]	sqrt_tbl	Precomputed array of square roots. The table element at index i should contain the value sqrt(i). Values from 0 to 2*el+1 must be precomputed (i.e. sqrt_tbl should contian 2*el+2 elements).

Return values

none

Author

Jason McEwen

ssht_dl_halfpi_trapani_fill_eighth2quarter_table()

void ssht_dl_halfpi_trapani_fill_eighth2quarter_table	(	double *	dl,
		int	L,
		ssht_dl_size_t	dl_size,
		int	el,
		double *	signs
	)

Fill in quarter Wigner plane for m = 0:l and mm = 0:l from the eighth m = 0:l and mm = 0:m. Takes a table of precomputed signs to avoid recomputing them.

Parameters

[in,out]	dl	Wigner plane. On input should contain values computed for one eighth of plane. On output will contain values computed for quarter of plane, using symmetries to fill in these values.
[in]	L	Harmonic band-limit.
[in]	dl_size	Size type of the memory to allocate.
[in]	el	Harmonic index to compute Wigner plane for.
[in]	signs	Precomputed array of signs. The array element at index i should contain the value (-1)^i. Values from 0 to el must be precomputed (i.e. signs should contian el+1 elements).

Return values

none

Author

Jason McEwen

ssht_dl_halfpi_trapani_fill_eighth2righthalf_table()

void ssht_dl_halfpi_trapani_fill_eighth2righthalf_table	(	double *	dl,
		int	L,
		ssht_dl_size_t	dl_size,
		int	el,
		double *	signs
	)

Fill in half Wigner plane for m = 0:l and mm = -l:l from the eighth m = 0:l and mm = 0:m. Takes a table of precomputed signs to avoid recomputing them.

Parameters

[in,out]	dl	Wigner plane. On input should contain values computed for one eighth of plane. On output will contain values computed for right half of plane, using symmetries to fill in these values.
[in]	L	Harmonic band-limit.
[in]	dl_size	Size type of the memory to allocate.
[in]	el	Harmonic index to compute Wigner plane for.
[in]	signs	Precomputed array of signs. The array element at index i should contain the value (-1)^i. Values from 0 to el must be precomputed (i.e. signs should contian el+1 elements).

Return values

none

Author

Jason McEwen

ssht_dl_halfpi_trapani_quarter_table()

void ssht_dl_halfpi_trapani_quarter_table	(	double *	dl,
		int	L,
		ssht_dl_size_t	dl_size,
		int	el,
		double *	sqrt_tbl
	)

Calculates quarter (for m = 0:l and mm = 0:l) of lth plane of a d-matrix for PI/2 using Trapani & Navaza's recursion method. For l>0, require the dl plane to be computed already with values for l-1. Also takes a table of precomputed square roots of integers to avoid recomputing them.

Warning

THIS GOES UNSTABLE SOMEWHERE ABOVE A BAND-LIMIT OF 1024. Best to compute eighth of plane and then extend to quarter using symmetry (see ssht_dl_halfpi_trapani_eighth_table and ssht_dl_halfpi_trapani_fill_eighth2quarter_table).

Parameters

[in,out]	dl	Wigner plane. On input this should be initialised to the plane computed for el-1. On output this will be replaced with the computed plane for el.
[in]	L	Harmonic band-limit.
[in]	dl_size	Size type of the memory to allocate.
[in]	el	Harmonic index to compute Wigner plane for.
[in]	sqrt_tbl	Precomputed array of square roots. The table element at index i should contain the value sqrt(i). Values from 0 to 2*el+1 must be precomputed (i.e. sqrt_tbl should contian 2*el+2 elements).

Return values

none

Author

Jason McEwen