primitive.cpp

普林斯顿开发的快速球面调和变换算法

// primitive.cpp: implementation of the primitive class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "primitive.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

primitive::primitive()
{

}

primitive::~primitive()
{

}

double primitive::L2_an(int m, int l)
{
	return (sqrt((((double) (2*l+3))/((double) (2*l+1))) *
	       (((double) (l-m+1))/((double) (l+m+1)))) *
	  (((double) (2*l+1))/((double) (l-m+1))));

}

double primitive::L2_cn(int m, int l)
{
	if (l != 0) {
    return (-1.0 *
	  sqrt((((double) (2*l+3))/((double) (2*l-1))) *
	       (((double) (l-m+1))/((double) (l+m+1))) *
	       (((double) (l-m))/((double) (l+m)))) *
	  (((double) (l+m))/((double) (l-m+1))));
  }
  else
    return 0.0;

}


double primitive::L2_cn_inv(int m, int l)
{
	  double dl, dm;

  dl = (double) l;
  dm = (double) m;

  return ( -(1.0 + (1. - 2. * dm)/(dm + dl)) *
	   sqrt( ((-1. + 2.*dl)/(3. + 2.*dl)) *
		 ((dl + dl*dl + dm + 2.*dl*dm + dm*dm)/
		  (dl + dl*dl - dm - 2.*dl*dm + dm*dm)) )
	   );


}



double primitive::L2_ancn(int m, int l)
{
	  double dl, dm;

  dl = (double) l;
  dm = (double) m;

  return( sqrt( 4.0 + ( (4.0 * dm * dm - 1.0)/
			(dl * dl - dm * dm) ) ) );

}

void primitive::vec_add(double *data1, double *data2, double *result, int n)
{
	  int k;


  for (k = 0; k < n % 4; ++k)
    result[k] = data1[k] + data2[k];

  for ( ; k < n ; k += 4)
    {
      result[k] = data1[k] + data2[k];
      result[k + 1] = data1[k + 1] + data2[k + 1];
      result[k + 2] = data1[k + 2] + data2[k + 2];
      result[k + 3] = data1[k + 3] + data2[k + 3];
    }

}


void primitive::vec_mul(double scalar, double *data1, double *result, int n)
{
	   int k;


   for( k = 0; k < n % 4; ++k)
     result[k] = scalar * data1[k];

   for( ; k < n; k +=4)
     {
       result[k] = scalar * data1[k];
       result[k + 1] = scalar * data1[k + 1];
       result[k + 2] = scalar * data1[k + 2];
       result[k + 3] = scalar * data1[k + 3];
     }


}

void primitive::vec_pt_mul(double *data1, double *data2, double *result, int n)
{
	   int k;

  
  for(k = 0; k < n % 4; ++k)
    result[k] = data1[k] * data2[k];
  
  for( ; k < n; k +=4)
    {
      result[k] = data1[k] * data2[k];
      result[k + 1] = data1[k + 1] * data2[k + 1];
      result[k + 2] = data1[k + 2] * data2[k + 2];
      result[k + 3] = data1[k + 3] * data2[k + 3];
    }

}

void primitive::ArcCosEvalPts(int n, double *eval_pts)
{
	    int i;
    double twoN;

    twoN = (double) (2 * n);

   for (i=0; i<n; i++)
     eval_pts[i] = (( 2.0*((double)i)+1.0 ) * PI) / twoN;

}


void primitive::EvalPts(int n, double *eval_pts)
{
	    int i;
    double twoN;

    twoN = (double) (2*n);

   for (i=0; i<n; i++)
     eval_pts[i] = cos((( 2.0*((double)i)+1.0 ) * PI) / twoN);


}



void primitive::Pmm_L2(int m, double *eval_pts, int n, double *result)
{
	  int i;
  double md, id, mcons;

  id = (double) 0.0;
  md = (double) m;
  mcons = sqrt(md + 0.5);

  for (i=0; i<m; i++) {
    mcons *= sqrt((md-(id/2.0))/(md-id));
    id += 1.0;
  }
  if (m != 0 )
    mcons *= pow(2.0,-md/2.0);
  if ((m % 2) != 0) mcons *= -1.0;

  for (i=0; i<n; i++) 
    result[i] = mcons * pow(sin(eval_pts[i]),((double) m));


}


void primitive::P_eval(int m, double *coeffs, double *eval_args, double *result, double *workspace, int bw)
{
	    double *prev, *prevprev, *temp1, *temp2, *temp3, *temp4, *x_i;
    int i, j, n;
    double splat;

    prevprev = workspace;
    prev = prevprev + (2*bw);
    temp1 = prev + (2*bw);
    temp2 = temp1 + (2*bw);
    temp3 = temp2 + (2*bw);
    temp4 = temp3 + (2*bw);
    x_i = temp4 + (2*bw);

    n = 2*bw;

    /* now get the evaluation nodes */
    this->EvalPts(n,x_i);

    /*   for(i=0;i<n;i++)
      fprintf(stderr,"in P_eval evalpts[%d] = %lf\n", i, x_i[i]);
      */   
    for (i=0; i<n; i++) 
      prevprev[i] = 0.0;

    if (m == 0) {
	for (i=0; i<n; i++) {
	  /* prev[i] = 0.707106781186547; sqrt(1/2) */
	   prev[i] = 1.0;
	    /* now mult by first coeff and add to result */
	    result[i] = coeffs[0] * prev[i];
	}
    }
    else {
	this->Pmm_L2(m, eval_args, n, prev);
	splat = coeffs[0];
	for (i=0; i<n; i++)
	  result[i] = splat * prev[i];
    }

    for (i=0; i<bw-m-1; i++) {
	this->vec_mul(L2_cn(m,m+i),prevprev,temp1,n);
	this->vec_pt_mul(prev, x_i, temp2, n);
	this->vec_mul(L2_an(m,m+i), temp2, temp3, n);
	this->vec_add(temp3, temp1, temp4, n); /* temp4 now contains P(m,m+i+1) */
	/* now add weighted P(m,m+i+1) to the result */
	splat = coeffs[i+1];
	for (j=0; j<n; j++)
	  result[j] += splat * temp4[j];
	memcpy(prevprev, prev, sizeof(double) * n);
	memcpy(prev, temp4, sizeof(double) * n);
    }

}