ccSHT is a C library comprised of routines that compute the discrete spherical harmonic transform (SHT). The SHT is essentially the Fourier transform for the surface of the unit sphere (S2). The basis functions in the SHT are an orthonormal set of eigenfunctions of the Laplace operator on S2 which span S2. The basis functions of the Fourier transform are an orthonormal set of eigenfunctions of the Laplace operator on n dimensional Euclidean space (Rn) which span Rn. The Fourier transform is ubiquitous and many libraries exist to do discrete fast Fourier transforms (FFT's) in Rn. There are many applications, however, where spherical coordinates arise naturally, and spectral analysis must be performed. In these cases the natural tool to use is the SHT.
A few other libraries exist which compute the discrete SHT but they apply to a more limited set of pixelizations and do not contain implementations for parallel architectures. ccSHT can be applied to a wide variety of pixelizations and contains both parallel and serial versions of the transform. The parallel code uses the MPI library to perform communications. ccSHT does the computation quickly (in O(N^(3/2)) operations) where N is the number of pixels.
There are a set of functions in the source code files ccSHTfortran.c and ccSHTmpiFortran.c which allow the primary functions of the ccSHT library to be easily linked from a FORTRAN code. For more information about the calling syntax of these FORTRAN compatible subroutines see the man page ccSHTfortran(1).
FFT's are used for part of the computation. These FFT's are done with the FFTW library which was chosen primarily for its portability. The Fourier transforms do not account for a significant amount of time in the computation (the time is proportional to the N * ln(N) where N is the number of pixels). None the less, FFTW is quite efficient on most architectures, so this short calculation is done quickly. For more information on the FFTW library visit their web site: http://www.fftw.org At the time of the release of ccSHT version 1.03, FFTW had just released version 3.0. This new version of FFTW has a different API, and so ccSHT must be used in conjunction with FFTW version 2.15. This code is now included with the ccSHT package in the form as downloaded from the fftw site on July 15, 2003.
The calculation of the SHT may be accelerated if the code is allowed to store some intermediate calculations on disk. These intermediate calculations will be stored if a directory called ccSHT_files exists in the working directory when a transform is performed. If the code can not find the file it is looking for, but the directory exists, then it will create the file. Using these files may speed up or slow down the calculation, depending on the disk access speed. In particular, when doing the computation on many processors, it may be faster not to use the files since many processors will need access to the disk at the same time.
ccSHT was designed to accommodate as many pixelizations as possible while still allowing the use of FFT's to speed the calculation. The information about the pixelization is stored in a structure called a coordStruct. This structure can be created from a list of right ascension and declination coordinates with the function convertRaDec2Coords(), and the structure is destroyed with the function called destroyCoords(). See the man pages of these two functions for more information.
In this documentation (and the ccSHT code) standard mathematical spherical coordinates (theta, phi) are used. We will briefly define these coordinates. Let us refer to an arbitrary fixed direction as the ``north pole''. theta is the latitudinal angle measured in radians with the value zero at the north pole and increasing to pi at the south pole. phi is the longitudinal angle measured in radians and increasing counter clockwise (to the east) around the north pole ranging from zero to two pi.
The pixelization must conform to the following list of constraints:
As a result of our loose restrictions on the pixelization of the sphere the data structure which describes the pixelization is somewhat complex. The best way to create this structure is to use convertRaDec2Coords(). A description of the structure is provided here for completeness, but it is not necessary to understand the details of the data structure to be able to use the transforms.
The below description assumes a particular ordering for the pixels. This ordering is used to vectorize discrete functions on the sphere. The ordering must conform to the following specification: two pixels which have the same theta value and are contiguous on the sphere (or only separated by a gap in the pixelization) must be contiguous in the pixel ordering unless they are the first and last pixels listed for a particular value of theta. This specification does not completely determine an ordering for a pixelization. The other choices which must be made to order the pixelization are arbitrary, but what ever ordering is chosen must be consistent throughout.
The structure coordStruct has the following fields:
[ 2 5 9 ]
[ 3 1 2 ]
and the input pixel vector was
[ 0.1 1 2 3 4 5 6 7 8 9 10 11 ]
then the zero padded vector (i.e. the vector with filled in gaps) would be
[ 0.1 1 0 0 0 2 3 4 0 5 6 7 8 0 0 9 10 11 ]
The spherical harmonic coefficients (a_{l,m}) and the scaled associated Legendre polynomials (Q_{l,m}) are indexed by l and m. For the spherical harmonic coefficients l ranges over the integers from zero to some band limit (lmax), and for each l, m ranges over the integers from -l to l. The associated Legendre polynomials are defined in a similar way for l from 0 to lmax, and m from 0 to l.
In order to store a set of spherical harmonic coefficients in memory we must define a way of vectorizing the set. To insure efficient memory access we have chosen to index in l major order since the calculation is done with this ordering. In the header file ccSHT.h there is a macro called almIndex which takes three integers (l,m,lmax) and returns a linear indexing. This ordering is l major, and increasing with both l and m and is defined below (note that lmax has been replaced by L for clarity):
almIndex(l,m,L) = (L*(L+1+2*m)+m*(2-abs(m+1)))/2+l
This implies that an a_{l,m} vector would be written:
[ a_{L,-L}, a_{L-1,-L+1}, a_{L,-L+1}, a_{L-2,-L+2}, ... , a_{0,0}, a_{1,0}, ... , a_{L,L} ]
This is the ordering for the output of forwardSHT and the input for backwardSHT.
There is a function for computing scaled associated Legendre polynomials called calculateQlm(). The data vector created by this function also uses an indexing scheme which is l major and increasing with both l and m. It is also defined in the header file ccSHT.h with the macro called Qindex which is defined:
Qindex(l,m,L) = L*m - (m*(m-1))/2 + l
Therefore a Q_{l,m} vector would be written:
[ Q_{0,0}, Q_{1,0}, ... , Q_{L,0}, Q_{1,1}, Q_{2,1}, ... , Q_{L, L} ]
This is the ordering for the output of calculateQlm().
Copyright (C) 2003 Christopher M. Cantalupo
ccSHT is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
Christopher Cantalupo <cmc@nersc.gov>
Send bug reports or comments to the above address.
forwardSHT(1), backwardSHT(1), forwardSHTmpi(1), backwardSHTmpi(1), calculateQlm(1), lmCalculations(1), convertRaDec2Coords(1), destroyCoords(1), fftw(1), ccSHTfortran(1)