/***************************************************************************
  **************************************************************************
  
                           S2kit 1.0

          A lite version of Spherical Harmonic Transform Kit

   Peter Kostelec, Dan Rockmore
   {geelong,rockmore}@cs.dartmouth.edu
  
   Contact: Peter Kostelec
            geelong@cs.dartmouth.edu
  
   Copyright 2004 Peter Kostelec, Dan Rockmore

   This file is part of S2kit.

   S2kit 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.

   S2kit 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 S2kit; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

   See the accompanying LICENSE file for details.
  
  ************************************************************************
  ************************************************************************/


/*************************************************************************/

/* Source code to synthesize functions using a naive method
   based on recurrence.  This is slow but does not require any
   precomputed functions, and is also stable. 
*/


#include <math.h>
#include <string.h>


/************************************************************************/
/*

   Naive_AnalysisX: computing the discrete Legendre transform of
                    a function via summing naively. I.e. This is
		    the FORWARD discrete Legendre transform.

   bw - bandwidth
   m - order

   data - a pointer to double array of size (2*bw) containing
          the sample points

   result - a pointer to double array of size (bw-m) which, at the
            conclusion of the routine, will contains the coefficients

   plmtable - a pointer to a double array of size (2*bw*(bw-m));
	      contains the PRECOMPUTED plms, i.e. associated Legendre
	      functions.  E.g. Should be generated by a call to

	      PmlTableGen()

	      (see pmls.c)

	      NOTE that these Legendres are normalized with norm
	      equal to 1 !!!

   workspace - array of size 2 * bw;


*/


void Naive_AnalysisX(double *data,
		     int bw,
		     int m,
		     double *weights,
		     double *result,
		     double *plmtable,
		     double *workspace)
{
  int i, j;
  double result0, result1, result2, result3;
  register double *wdata;

  wdata = workspace;

  /* make sure result is zeroed out */
  memset( result, 0, sizeof(double) * (bw - m) );

  /* apply quadrature weights */
  /*
    I only have to differentiate between even and odd
    weights when doing something like seminaive, something
    which involves the dct. In this naive case, the parity of
    the order of the transform doesn't matter because I'm not
    dividing by sin(x) when precomputing the Legendres (because
    I'm not taking their dct). The plain ol' weights are just
    fine. 
  */
  
  for(i = 0; i < 2 * bw; i++)
    wdata[i] = data[i] * weights[i];

  /* unrolling seems to work */
  if ( 1 )
    {
      for (i = 0; i < bw - m; i++)
	{
	  result0 = 0.0; result1 = 0.0;
	  result2 = 0.0; result3 = 0.0;
	  
	  for(j = 0; j < (2 * bw) % 4; ++j)
	    result0 += wdata[j] * plmtable[j];

	  for( ; j < (2 * bw); j += 4)
	    {
	      result0 += wdata[j] * plmtable[j];
	      result1 += wdata[j + 1] * plmtable[j + 1];
	      result2 += wdata[j + 2] * plmtable[j + 2];
	      result3 += wdata[j + 3] * plmtable[j + 3];
	    }
	  result[i] = result0 + result1 + result2 + result3;

	  plmtable += (2 * bw);
	}
    }
  else
    {
      for (i = 0; i < bw - m; i++)
	{
	  result0 = 0.0 ; 	  
	  for(j = 0; j < (2 * bw) ; j++)
	    result0 += wdata[j] * plmtable[j];
	  result[i] = result0 ;

	  plmtable += (2 * bw);
	}
    }
}


/************************************************************************/
/* This is the procedure that synthesizes a function from a list
   of coefficients of a Legendre series. I.e. this is the INVERSE
   discrete Legendre transform.

   Function is synthesized at the (2*bw) Chebyshev nodes. Associated
   Legendre functions are assumed to be precomputed.
   
   bw - bandwidth

   m - order

   coeffs - a pointer to double array of size (bw-m).  First coefficient is
            coefficient for Pmm
   result - a pointer to double array of size (2*bw); at the conclusion
            of the routine, this array will contain the
            synthesized function

   plmtable - a pointer to a double array of size (2*bw*(bw-m));
	      contains the PRECOMPUTED plms, i.e. associated Legendre
	      functions. E.g. Should be generated by a call to

	      PmlTableGen(),

	      (see pmls.c)


	      NOTE that these Legendres are normalized with norm
	      equal to 1 !!!
*/

void Naive_SynthesizeX(double *coeffs,
		       int bw,
		       int m,
		       double *result,
		       double *plmtable)
{
  int i, j;
  double tmpcoef;

  /* make sure result is zeroed out */
  memset( result, 0, sizeof(double) * 2 * bw );

  for ( i = 0 ; i < bw - m ; i ++ )
    {
      tmpcoef = coeffs[i] ;
      if ( tmpcoef != 0.0 )
	for (j=0; j<(2*bw); j++)
	  result[j] += (tmpcoef * plmtable[j]);
      plmtable += (2 * bw ) ;
    }
}
