/***************************************************************************  **************************************************************************                  Spherical Harmonic Transform Kit 2.7       Contact: Peter Kostelec            geelong@cs.dartmouth.edu       Copyright 1997-2003  Sean Moore, Dennis Healy,                        Dan Rockmore, Peter Kostelec     Copyright 2004  Peter Kostelec, Dan Rockmore         SpharmonicKit 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.       SpharmonicKit 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.       You should have received a copy of the GNU General Public License     along with this program; if not, write to the Free Software     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.       Commercial use is absolutely prohibited.     See the accompanying LICENSE file for details.    ************************************************************************  ************************************************************************/==================================================February, 2004: SpharmonicKit 2.7 releasedAnother (and hopefully last) scaling discrepancy corrected.A user has identified another scaling error, this timein the TransMult() routine which is used for convolutions.The coefficients were not being weighted correctly prior toperforming the inverse spherical transform.In the original Driscoll-Healy paper, in Theorem 1 it is provedthat, for two function f and h defined on the sphere, thetransform of their convolution is a pointwise product ofthe 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 thecode.We thank the user for identifying this error, and providing thefollowing 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 releasedScaling discrepancy corrected.While the code is internally consistent, as far as normalizations areconcerned, a user recently pointed out to us that our scaling of theY_lm's differs from the usual Y_lm's.If, for example, one sampled the spherical harmonic Y_1^0 on thesphere in Mathematica (tm), and then used SpharmonicKit to takeits forward spherical transform, the computed coefficientwould not be 1, but rather 1/(2*sqrt(PI)). And if one sampledY_1^1, the computed coefficient would 1/sqrt(2*PI), and not 1.So, depending on the order of the spherical transform, we were missingeither a 2*sqrt(PI), or a sqrt(2*PI). We are grateful for the userpointing this out to us, and apologize for any inconvenience this mayhave caused. The spherical transform routines, including those havingto do with convolution, have been corrected in version 2.6.To detail the fix:A) Forward Spherical TransformTaking 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 computedcoefficients are in the arrays "rresult" and "iresult", do somethinglike 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 TransformTaking 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 coefficientsare in the arrays "rcoeffs" and "icoeffs", then you should do somethinglike 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 makeeverything a (quicker) multiply.=========================================July, 1998: SpharmonicKit 2.5 releasedThis version of the Kit is designed to use a slight variationof FFTPACK, the freely available collection of FORTRAN subprogramsfor "... calculating fast Fourier transforms for both complex andreal periodic sequences and certain other symmetric sequences ...."To be precise, the Kit now uses slight modifications of the FFTand DCT routines found in FFTPACK. These routines are considerablymore efficient than those provided in SpharmonicKit 2 (and are stillincluded in this current release). Details can be found in HOWTO_FFTPACK.Making the reasonable assumption that the input data will alwaysbe strictly real, the fft-portions of all the spherical transformswere modified to take full advantage of the symmetries that thisassumption brings. This is only in the case when FFTPACK is used.Up to now, the symmetries were only partially and not fullyexploited.Also redesigned some of the documentation.==========================================March, 1998: SpharmonicKit 2 releasedIn the original release of the Kit, all the Legendre transforms(and code relying on Legendre transforms) were based on theseminaive and naive algorithms. Code based on the work of Driscolland 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 arethe 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 1024The seminaive spherical algorithms (forward and reverse)were modified and versions which read precomputed dataoff disk are provided, as well.