ccSHT(1) ccSHT(1) NAME ccSHT - A Spherical Harmonic Transform Library DESCRIPTION Introduction 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 opera- tor 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. Other SHT Libraries 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 per- form communications. ccSHT does the computation quickly (in O(N^(3/2)) operations) where N is the number of pix- els. Using ccSHT with FORTRAN There are a set of functions in the source code files ccSHTfortran.c and ccSHTmpiFortran.c which allow the pri- mary functions of the ccSHT library to be easily linked from a FORTRAN code. For more information about the call- ing syntax of these FORTRAN compatible subroutines see the man page ccSHTfortran(1). Fourier Transforms 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 computa- tion (the time is proportional to the N * ln(N) where N is the number of pixels). None the less, FFTW is quite effi- cient 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. Intermediate Calculation Files 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 direc- tory 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. Pixelization 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 conver- tRaDec2Coords(), 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 mathe- matical 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: + The pixels must lay on rows of constant latitude (theta). + Quadrature weights must be a constant. If the quadrature weights are not constant, but a function only of the pixel location, then the input vector should be point wise multiplied by by the corre- sponding quadrature weight vector before a forward transform is computed. + The longitudinal (phi) separation between pixels must be constant for each row of fixed theta, but different rows can have different longitudinal pixel separations. + There may be missing pixels (called gaps) within rows of constant latitude (theta), but the underly- ing pixel spacing must be even. Another way of saying this is that entire pixels may be unspeci- fied but not fractions of pixels. + The pixelization is not required to cover the entire sphere. In particular, rows of pixels are not required to entirely encircle the sphere. Every row of pixels has a first and last pixel. These pixels may or may not be contiguous on the sphere. If they are not contiguous then there is a space before the first pixel and after the last pixel where no pixels are specified. This space is not what we refer to as a gap which refers to miss- ing pixels in a row after the first pixel and before the last pixel. The region before the first pixel and after the last pixel in a row is in effect a special gap. The function value at all unspecified regions (gaps) in the pixelization is inferred to be zero. The coordStruct Data Structure As a result of our loose restrictions on the pixelization of the sphere the data structure which describes the pix- elization is somewhat complex. The best way to create this structure is to use convertRaDec2Coords(). A description of the structure is provided here for com- pleteness, but it is not necessary to understand the details of the data structure to be able to use the trans- forms. 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 sep- arated by a gap in the pixelization) must be contiguous in the pixel ordering unless they are the first and last pix- els listed for a particular value of theta. This specifi- cation 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 order- ing is chosen must be consistent throughout. The structure coordStruct has the following fields: nPix Number of specified pixels on the sphere. nThetaVals Number of latitudinal (constant theta) rows of pix- els which are specified. thetaVals A vector of theta values in radians for each row of pixels. This list should be ordered in a way con- sistent with the pixel indexing scheme. thetaBreaks A vector of pixel numbers for the start of each row of constant theta. phi0 A vector of phi values for the first pixel in each latitudinal (constant theta) row. deltaPhi A vector of pixel widths in radians for each row of constant theta. This is the change in phi between two contiguous pixels with the same theta value. The sign of each deltaPhi specifies if the phi val- ues increase or decrease (+/-) in the pixel order- ing for the corresponding row. Note that deltaPhi is not the angular distance between contiguous pix- els which is deltaPhi*sin(theta). nGaps Number of gaps in otherwise regularly spaced in phi rows of constant theta. Note that these gaps must have a width which is an integer number of pixels, and if there is an unspecified region before the first pixel in a row and after the last pixel in the row, this is not considered a gap. gaps A column major ordered 2 by nGaps matrix of pixel numbers where gaps begin in the first row and width of gap in second row. So if gaps were the follow- ing matrix [ 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 ] pixSize Quadrature weight for use in the forward transform (i.e. area of pixels in radians^2). If quadrature weight is not a constant, set pixSize to one (in which case the integral in the forward SHT is just a sum) and point wise multiply your input vector by the desired quadrature weighting scheme before for- ward transforming. Data Structures for the Spherical Harmonic Coefficients and Legendre Polynomials 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 coeffi- cients 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 Legen- dre 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 Version 1.03 July 2003 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 Send bug reports or comments to the above address. SEE ALSO forwardSHT(1), backwardSHT(1), forwardSHTmpi(1), backward- SHTmpi(1), calculateQlm(1), lmCalculations(1), conver- tRaDec2Coords(1), destroyCoords(1), fftw(1), ccSHTfor- tran(1) Version 1.03: July 2003 ccSHT(1)