trunk/SHAPES/primitive.c

/***************************************************************************
  **************************************************************************
  
                           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.
  
  ************************************************************************
  ************************************************************************/


/*

  some "primitive" functions that are used in cospmls.c
  and newmathx.c

  */

#include <math.h>
#include <string.h>  /* to declare memcpy */

#ifndef PI
#define PI 3.14159265358979
#endif

/************************************************************************/
/* Recurrence coefficients */
/************************************************************************/
/* Recurrence coefficents for L2-normed associated Legendre
   recurrence.  When using these coeffs, make sure that
   inital Pmm function is also L2-normed */
/* l represents degree, m is the order */

double 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))));

}

/* note - if input l is zero, need to return 0 */
double 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;

}

/* when using the reverse recurrence, instead of calling
   1/L2_cn_tr(m,l), let me just define the function ...
   it might be more stable */

double 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)) )
           );

}

/* when using the reverse recurrence, instead of calling
   -L2_an(m,l)/L2_cn(m,l), let me just define the
   function ... it might be more stable */

double 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) ) ) );
}

/************************************************************************/
/* vector arithmetic operations */
/************************************************************************/
/* does result = data1 + data2 */
/* result and data are vectors of length n */

void 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];
    }
}
/************************************************************************/
/************************************************************************/
/*
   vec_mul(scalar,data1,result,n) multiplies the vector 'data1' by
   'scalar' and returns in result 
*/
void 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];
     }

}
/************************************************************************/
/* point-by-point multiplication of vectors */

void 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];
    }
 
}


/************************************************************************/
/* returns an array of the angular arguments of n Chebyshev nodes */
/* eval_pts points to a double array of length n */

void 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;

}
/************************************************************************/
/* returns an array of n Chebyshev nodes */

void 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);

}

/************************************************************************/
/* L2 normed Pmm.  Expects input to be the order m, an array of
 evaluation points arguments of length n, and a result vector of length n */
/* The norming constant can be found in Sean's PhD thesis */
/* This has been tested and stably computes Pmm functions thru bw=512 */

void 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));

}

/************************************************************************/
/************************************************************************/
/*
  This piece of code synthesizes a function which is the weighted sum of 
  L2-normalized associated Legendre functions.

  The coeffs array should contain bw - m coefficients ordered from
  zeroth degree to bw-1, and eval_pts should be an array of the
  arguments (arccos) of the desired evaluation points of length 2*bw. 

  Answer placed in result (and has length 2*bw).
  
  workspace needs to be of size 16 * bw
  
*/

void 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 */
    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) */

          /* now mult by first coeff and add to result */
          result[i] = coeffs[0] * prev[i];
        }
    }
    else {
        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++) {
        vec_mul(L2_cn(m,m+i),prevprev,temp1,n);
        vec_pt_mul(prev, x_i, temp2, n);
        vec_mul(L2_an(m,m+i), temp2, temp3, n);
        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);
    }

}

Revision:	1287
Committed:	Wed Jun 23 20:18:48 2004 UTC (20 years ago) by chrisfen
Content type:	text/plain
File size:	8525 byte(s)
Log Message:	Major progress towards inclusion of spherical harmonic transform capability - still having some build issues...
#	User	Rev	Content
1	chrisfen	1287	/***************************************************************************
2			**************************************************************************
3
4			S2kit 1.0
5
6			A lite version of Spherical Harmonic Transform Kit
7
8			Peter Kostelec, Dan Rockmore
9			{geelong,rockmore}@cs.dartmouth.edu
10
11			Contact: Peter Kostelec
12			geelong@cs.dartmouth.edu
13
14			Copyright 2004 Peter Kostelec, Dan Rockmore
15
16			This file is part of S2kit.
17
18			S2kit is free software; you can redistribute it and/or modify
19			it under the terms of the GNU General Public License as published by
20			the Free Software Foundation; either version 2 of the License, or
21			(at your option) any later version.
22
23			S2kit is distributed in the hope that it will be useful,
24			but WITHOUT ANY WARRANTY; without even the implied warranty of
25			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26			GNU General Public License for more details.
27
28			You should have received a copy of the GNU General Public License
29			along with S2kit; if not, write to the Free Software
30			Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
31
32			See the accompanying LICENSE file for details.
33
34			************************************************************************
35			************************************************************************/
36
37
38
39			/*
40
41			some "primitive" functions that are used in cospmls.c
42			and newmathx.c
43
44			*/
45
46			#include <math.h>
47			#include <string.h> /* to declare memcpy */
48
49			#ifndef PI
50			#define PI 3.14159265358979
51			#endif
52
53			/************************************************************************/
54			/* Recurrence coefficients */
55			/************************************************************************/
56			/* Recurrence coefficents for L2-normed associated Legendre
57			recurrence. When using these coeffs, make sure that
58			inital Pmm function is also L2-normed */
59			/* l represents degree, m is the order */
60
61			double L2_an(int m,
62			int l)
63			{
64			return (sqrt((((double) (2l+3))/((double) (2l+1))) *
65			(((double) (l-m+1))/((double) (l+m+1)))) *
66			(((double) (2*l+1))/((double) (l-m+1))));
67
68			}
69
70			/* note - if input l is zero, need to return 0 */
71			double L2_cn(int m,
72			int l)
73			{
74			if (l != 0) {
75			return (-1.0 *
76			sqrt((((double) (2l+3))/((double) (2l-1))) *
77			(((double) (l-m+1))/((double) (l+m+1))) *
78			(((double) (l-m))/((double) (l+m)))) *
79			(((double) (l+m))/((double) (l-m+1))));
80			}
81			else
82			return 0.0;
83
84			}
85
86			/* when using the reverse recurrence, instead of calling
87			1/L2_cn_tr(m,l), let me just define the function ...
88			it might be more stable */
89
90			double L2_cn_inv(int m,
91			int l)
92			{
93			double dl, dm;
94
95			dl = (double) l;
96			dm = (double) m;
97
98			return ( -(1.0 + (1. - 2. * dm)/(dm + dl)) *
99			sqrt( ((-1. + 2.dl)/(3. + 2.dl)) *
100			((dl + dldl + dm + 2.dldm + dmdm)/
101			(dl + dldl - dm - 2.dldm + dmdm)) )
102			);
103
104			}
105
106			/* when using the reverse recurrence, instead of calling
107			-L2_an(m,l)/L2_cn(m,l), let me just define the
108			function ... it might be more stable */
109
110			double L2_ancn(int m,
111			int l)
112			{
113			double dl, dm;
114
115			dl = (double) l;
116			dm = (double) m;
117
118			return( sqrt( 4.0 + ( (4.0 * dm * dm - 1.0)/
119			(dl * dl - dm * dm) ) ) );
120			}
121
122			/************************************************************************/
123			/* vector arithmetic operations */
124			/************************************************************************/
125			/* does result = data1 + data2 */
126			/* result and data are vectors of length n */
127
128			void vec_add(double *data1,
129			double *data2,
130			double *result,
131			int n)
132			{
133			int k;
134
135
136			for (k = 0; k < n % 4; ++k)
137			result[k] = data1[k] + data2[k];
138
139			for ( ; k < n ; k += 4)
140			{
141			result[k] = data1[k] + data2[k];
142			result[k + 1] = data1[k + 1] + data2[k + 1];
143			result[k + 2] = data1[k + 2] + data2[k + 2];
144			result[k + 3] = data1[k + 3] + data2[k + 3];
145			}
146			}
147			/************************************************************************/
148			/************************************************************************/
149			/*
150			vec_mul(scalar,data1,result,n) multiplies the vector 'data1' by
151			'scalar' and returns in result
152			*/
153			void vec_mul(double scalar,
154			double *data1,
155			double *result,
156			int n)
157			{
158			int k;
159
160
161			for( k = 0; k < n % 4; ++k)
162			result[k] = scalar * data1[k];
163
164			for( ; k < n; k +=4)
165			{
166			result[k] = scalar * data1[k];
167			result[k + 1] = scalar * data1[k + 1];
168			result[k + 2] = scalar * data1[k + 2];
169			result[k + 3] = scalar * data1[k + 3];
170			}
171
172			}
173			/************************************************************************/
174			/* point-by-point multiplication of vectors */
175
176			void vec_pt_mul(double *data1,
177			double *data2,
178			double *result,
179			int n)
180			{
181			int k;
182
183
184			for(k = 0; k < n % 4; ++k)
185			result[k] = data1[k] * data2[k];
186
187			for( ; k < n; k +=4)
188			{
189			result[k] = data1[k] * data2[k];
190			result[k + 1] = data1[k + 1] * data2[k + 1];
191			result[k + 2] = data1[k + 2] * data2[k + 2];
192			result[k + 3] = data1[k + 3] * data2[k + 3];
193			}
194
195			}
196
197
198			/************************************************************************/
199			/* returns an array of the angular arguments of n Chebyshev nodes */
200			/* eval_pts points to a double array of length n */
201
202			void ArcCosEvalPts(int n,
203			double *eval_pts)
204			{
205			int i;
206			double twoN;
207
208			twoN = (double) (2 * n);
209
210			for (i=0; i<n; i++)
211			eval_pts[i] = (( 2.0((double)i)+1.0 ) PI) / twoN;
212
213			}
214			/************************************************************************/
215			/* returns an array of n Chebyshev nodes */
216
217			void EvalPts( int n,
218			double *eval_pts)
219			{
220			int i;
221			double twoN;
222
223			twoN = (double) (2*n);
224
225			for (i=0; i<n; i++)
226			eval_pts[i] = cos((( 2.0((double)i)+1.0 ) PI) / twoN);
227
228			}
229
230			/************************************************************************/
231			/* L2 normed Pmm. Expects input to be the order m, an array of
232			evaluation points arguments of length n, and a result vector of length n */
233			/* The norming constant can be found in Sean's PhD thesis */
234			/* This has been tested and stably computes Pmm functions thru bw=512 */
235
236			void Pmm_L2( int m,
237			double *eval_pts,
238			int n,
239			double *result)
240			{
241			int i;
242			double md, id, mcons;
243
244			id = (double) 0.0;
245			md = (double) m;
246			mcons = sqrt(md + 0.5);
247
248			for (i=0; i<m; i++) {
249			mcons *= sqrt((md-(id/2.0))/(md-id));
250			id += 1.0;
251			}
252			if (m != 0 )
253			mcons *= pow(2.0,-md/2.0);
254			if ((m % 2) != 0) mcons *= -1.0;
255
256			for (i=0; i<n; i++)
257			result[i] = mcons * pow(sin(eval_pts[i]),((double) m));
258
259			}
260
261			/************************************************************************/
262			/************************************************************************/
263			/*
264			This piece of code synthesizes a function which is the weighted sum of
265			L2-normalized associated Legendre functions.
266
267			The coeffs array should contain bw - m coefficients ordered from
268			zeroth degree to bw-1, and eval_pts should be an array of the
269			arguments (arccos) of the desired evaluation points of length 2*bw.
270
271			Answer placed in result (and has length 2*bw).
272
273			workspace needs to be of size 16 * bw
274
275			*/
276
277			void P_eval(int m,
278			double *coeffs,
279			double *eval_args,
280			double *result,
281			double *workspace,
282			int bw)
283			{
284			double prev, prevprev, temp1, temp2, temp3, temp4, *x_i;
285			int i, j, n;
286			double splat;
287
288			prevprev = workspace;
289			prev = prevprev + (2*bw);
290			temp1 = prev + (2*bw);
291			temp2 = temp1 + (2*bw);
292			temp3 = temp2 + (2*bw);
293			temp4 = temp3 + (2*bw);
294			x_i = temp4 + (2*bw);
295
296			n = 2*bw;
297
298			/* now get the evaluation nodes */
299			EvalPts(n,x_i);
300
301			/* for(i=0;i<n;i++)
302			fprintf(stderr,"in P_eval evalpts[%d] = %lf\n", i, x_i[i]);
303			*/
304			for (i=0; i<n; i++)
305			prevprev[i] = 0.0;
306
307			if (m == 0) {
308			for (i=0; i<n; i++) {
309			prev[i] = 0.707106781186547; /* sqrt(1/2) */
310
311			/* now mult by first coeff and add to result */
312			result[i] = coeffs[0] * prev[i];
313			}
314			}
315			else {
316			Pmm_L2(m, eval_args, n, prev);
317			splat = coeffs[0];
318			for (i=0; i<n; i++)
319			result[i] = splat * prev[i];
320			}
321
322			for (i=0; i<bw-m-1; i++) {
323			vec_mul(L2_cn(m,m+i),prevprev,temp1,n);
324			vec_pt_mul(prev, x_i, temp2, n);
325			vec_mul(L2_an(m,m+i), temp2, temp3, n);
326			vec_add(temp3, temp1, temp4, n); /* temp4 now contains P(m,m+i+1) */
327			/* now add weighted P(m,m+i+1) to the result */
328			splat = coeffs[i+1];
329			for (j=0; j<n; j++)
330			result[j] += splat * temp4[j];
331			memcpy(prevprev, prev, sizeof(double) * n);
332			memcpy(prev, temp4, sizeof(double) * n);
333			}
334
335			}
336