================================================== February, 2004: SpharmonicKit 2.7 released Another (and hopefully last) scaling discrepancy corrected. A user has identified another scaling error, this time in the TransMult() routine which is used for convolutions. The coefficients were not being weighted correctly prior to performing the inverse spherical transform. In the original Driscoll-Healy paper, in Theorem 1 it is proved that, for two function f and h defined on the sphere, the transform of their convolution is a pointwise product of the transforms (f*h)^(l,m) = 2*pi*sqrt(4*pi/(2*l+1)) f^(l,m) * h^(l,0) where ^ denotes the Fourier coefficient of degree l, order m. In a nutshell, we forgot the 2*pi*sqrt(4*pi/(2*l+1)) in the code. We thank the user for identifying this error, and providing the following correction. The main for-loops in TransMult() used to be ***************************************** for (m=0; m‹bw; m++) { for (l=m; l‹bw; l++) { compmult(rfiltercoeffs[l], ifiltercoeffs[l], rdptr[l-m], idptr[l-m], rrptr[l-m], irptr[l-m]); } rdptr += bw-m; idptr += bw-m; rrptr += bw-m; irptr += bw-m; } for (m=bw+1; m‹size; m++) { for (l=size-m; l‹bw; l++){ compmult(rfiltercoeffs[l], ifiltercoeffs[l], rdptr[l-size+m], idptr[l-size+m], rrptr[l-size+m], irptr[l-size+m]); } rdptr += m-bw; idptr += m-bw; rrptr += m-bw; irptr += m-bw; } ***************************************** The corrected version is ***************************************** for (m=0; m‹bw; m++) { for (l=m; l‹bw; l++) { compmult(rfiltercoeffs[l], ifiltercoeffs[l], rdptr[l-m], idptr[l-m], rrptr[l-m], irptr[l-m]); rrptr[l-m] *= sqrt(4*M_PI/(2*l+1)); irptr[l-m] *= sqrt(4*M_PI/(2*l+1)); } rdptr += bw-m; idptr += bw-m; rrptr += bw-m; irptr += bw-m; } for (m=bw+1; m‹size; m++) { for (l=size-m; l‹bw; l++){ compmult(rfiltercoeffs[l], ifiltercoeffs[l], rdptr[l-size+m], idptr[l-size+m], rrptr[l-size+m], irptr[l-size+m]); rrptr[l-size+m] *= sqrt(4*M_PI/(2*l+1)); irptr[l-size+m] *= sqrt(4*M_PI/(2*l+1)); } rdptr += m-bw; idptr += m-bw; rrptr += m-bw; irptr += m-bw; } ***************************************** ================================================== July, 2003: SpharmonicKit 2.6 released Scaling discrepancy corrected. While the code is internally consistent, as far as normalizations are concerned, a user recently pointed out to us that our scaling of the Y_lm's differs from the usual Y_lm's. If, for example, one sampled the spherical harmonic Y_1^0 on the sphere in Mathematica (tm), and then used SpharmonicKit to take its forward spherical transform, the computed coefficient would not be 1, but rather 1/(2*sqrt(PI)). And if one sampled Y_1^1, the computed coefficient would 1/sqrt(2*PI), and not 1. Depending on the order of the spherical transform, we were missing either a 2*sqrt(PI), or a sqrt(2*PI). We are grateful for the user pointing this out to us, and apologize for any inconvenience this may have caused. The spherical transform routines, including those having to do with convolution, have been corrected in version 2.6. To detail the fix: A) Forward Spherical Transform Taking the forward spherical transform, as defined in version 2.5, the computed coefficients need to be scaled as follows: - multiply all order m = 0 coefficients by 2*sqrt(PI) - multiply all order m != 0 coefficients by sqrt(2*pi) For example, if the real and imaginary parts of the computed coefficients are in the arrays "rresult" and "iresult", do something like this: for(m=0;m‹bw;m++) for(l=m;l‹bw;l++){ dummy = seanindex(m,l,bw); if ( m == 0 ) { rresult[dummy] *= (2.*sqrt(PI)); iresult[dummy] *= (2.*sqrt(PI)); } else { rresult[dummy] *= (sqrt(2.*PI)); iresult[dummy] *= (sqrt(2.*PI)); } /* now for the negative-order coefficients */ if ( m != 0 ) { dummy = seanindex(-m,l,bw); rresult[dummy] *= (sqrt(2.*PI)); iresult[dummy] *= (sqrt(2.*PI)); } } B) Inverse Spherical Transform Taking the inverse spherical transform, as defined in version 2.5, the computed coefficients need to be scaled as follows, *BEFORE* applying the inverse spherical transform: - divide all order m = 0 coefficients by 2*sqrt(PI) - divide all order m != 0 coefficients by sqrt(2*pi) For example, if the real and imaginary parts of the computed coefficients are in the arrays "rcoeffs" and "icoeffs", then you should do something like this: for(m=0;m‹bw;m++) for(l=m;l‹bw;l++){ dummy = seanindex(m,l,bw); if ( m == 0 ) { rcoeffs[dummy] /= (2.*sqrt(PI)); icoeffs[dummy] /= (2.*sqrt(PI)); } else { rcoeffs[dummy] /= (sqrt(2.*PI)); icoeffs[dummy] /= (sqrt(2.*PI)); } /* now for the negative-order coefficients */ if ( m != 0 ) { dummy = seanindex(-m,l,bw); rcoeffs[dummy] /= (sqrt(2.*PI)); icoeffs[dummy] /= (sqrt(2.*PI)); } } To save some system calls to "sqrt", replace sqrt(2*PI), 2*sqrt(PI), 1/sqrt(2*PI) and 1/2*sqrt(PI) by their numerical equivalents and make everything a (quicker) multiply. ========================================= July, 1998: SpharmonicKit 2.5 released This version of the Kit is designed to use a slight variation of FFTPACK, the freely available collection of FORTRAN subprograms for "... calculating fast Fourier transforms for both complex and real periodic sequences and certain other symmetric sequences ...." To be precise, the Kit now uses slight modifications of the FFT and DCT routines found in FFTPACK. These routines are considerably more efficient than those provided in SpharmonicKit 2 (and are still included in this current release). Details can be found in HOWTO_FFTPACK. Making the reasonable assumption that the input data will always be strictly real, the fft-portions of all the spherical transforms were modified to take full advantage of the symmetries that this assumption brings. This is only in the case when FFTPACK is used. Up to now, the symmetries were only partially and not fully exploited. Also redesigned some of the documentation. ========================================== March, 1998: SpharmonicKit 2 released In the original release of the Kit, all the Legendre transforms (and code relying on Legendre transforms) were based on the seminaive and naive algorithms. Code based on the work of Driscoll and Healy is introduced in this release. As an extremely brief "what's new in this release" ... In this latest edition of the Kit, the following are the new major additions: Forward Legendre transforms: 1) a slight variation of the basic Driscoll-Healy (DH) Legendre transform algorithm; for bw = 16 through 1024 (must be power of 2) 2) the bounded DH-Mid algorithm; for bw = 16 through 1024 (must be power of 2) 3) the simple-split and hybrid algorithms; for bw = 16 through 1024 (must be power of 2) Forward Spherical Transforms: 1) a hybrid spherical transform (based on the hybrid Legendre transform) that precomputes in memory all necessary Legendre polynomial cosine transforms prior to transforming; for bw = 64 through 512 2) a hybrid spherical transform that reads the precomputed data off disk; for bw = 16 through 1024 3) a spherical convolution routine which uses the hybrid spherical transform in the forward direction and seminaive algorithm in the reverse; precomputes in memory prior to transforming; for bw = 64 through 512 4) a spherical convolution routine which uses the hybrid spherical transform in the forward direction and seminaive algorithm in the reverse; reads the precomputed data off disk; for bw = 16 through 1024 The seminaive spherical algorithms (forward and reverse) were modified and versions which read precomputed data off disk are provided, as well.

Return to README

Return to SpharmonicKit