 |
| |
|
| |
|
| |
|
| |
|
|
Asher
|
 |
Calgary, Alberta
Nov 1999
time: 22:37
|
|
quote:
Originally posted by Flip McWho
Ok so what does the OP actually say in layman? |
Isn't it obvious?
code:
/ **************************************************
**********************/
/***
In the relation
P_{L+r} = A_r^L P_L + B_r^L P_{L-1}
this routine generates B_r^L . Initial conditions are
B_0^L = 0, B_{-1}^L = 1.
m = order of the problem
r = how much to shift
l = where to start shift (i.e. l = L in above *)
n = number of points desired
result is double array of size n
workspace is double array of size 4 * n
***/
static void Bshifted( int m,
int r,
int l,
int n,
double *result,
double *workspace )
{
double *cospts, *pi, *piminus1, *pcurr;
double di, dn, aval, cval;
int i, j;
/* allocate work space */
cospts = workspace;
pi = cospts + n;
piminus1 = pi + n;
pcurr = piminus1 + n;
/* generate sample points */
dn = (double) n;
for (i=0; i<n; i++) {
di = (double) i;
cospts[i] = cos(
(((2.0 * di)+1.0) * PI)/(2.0 * dn));
}
if (r < 0) {
for (i=0; i<n; i++)
result[i] = 1.0;
}
/* piminus1 starts off as all zeros */
for (i=0; i<n; i++)
piminus1[i] = 0.0;
if (r == 0) {
for (i=0; i<n; i++)
result[i] = 0.0;
}
/* return (piminus1); */
/* generate initial value of pi */
for (i=0; i<n; i++)
pi[i] = L2_cn(m,m+l);
if (r == 1) {
for (i=0; i<n; i++)
result[i] = pi[i];
}
/* return(pi); */
/* main loop */
if (r > 1) {
for (i=m+1; i<m+r; i++) {
aval = L2_an(m,i+l);
cval = L2_cn(m,i+l);
for (j=0; j<n; j++)
pcurr[j] = aval * cospts[j] * pi[j] +
cval * piminus1[j];
/* need to copy values */
memcpy(piminus1, pi, sizeof(double) * n);
memcpy(pi, pcurr, sizeof(double) * n);
}
memcpy(result, pcurr, sizeof(double) * n);
}
/* return (pcurr); */
}
static void A_revshifted( int m,
int r,
int l,
int n,
double *result,
double *workspace )
{
double *cospts, *pi, *piplus1, *pcurr;
double di, dn, aval, cval;
int i, j;
/* allocate work space */
cospts = workspace;
pi = cospts + n;
piplus1 = pi + n;
pcurr = piplus1 + n;
/* generate sample points */
dn = (double) n;
for (i=0; i<n; i++) {
di = (double) i;
cospts[i] = cos((((2.0 * di)+1.0) * PI)/(2.0 * dn));
}
if (r < 0) {
for (i=0; i<n; i++)
result[i] = 0.0;
}
/* return (zeros(n)); */
/* piplus1 starts off as all ones */
for (i=0; i<n; i++)
piplus1[i] = 1.0;
if (r == 0) {
for (i=0; i<n; i++)
result[i] = 1.0;
}
/* return (piplus1); */
/* generate initial value of pi */
/* aval = -L2_an(m,m+l); */
/* cval = 1.0 / L2_cn(m,m+l); */
aval = L2_ancn(m,m+l);
for (i=0; i<n; i++)
pi[i] = cospts[i] * aval;
if (r == 1) {
for (i=0; i<n; i++)
result[i] = pi[i];
}
/* return(pi); */
/* main loop */
if (r > 1) {
for(i = 1 ; i < r ; i ++){
/* aval = -L2_an(m, m + l - i); */
aval = L2_ancn(m, m + l - i);
/* cval = 1.0/L2_cn(m, m + l - i); */
cval = L2_cn_inv(m,m+l-i);
for(j = 0 ; j < n ; j ++)
pcurr[j] = aval * cospts[j] * pi[j] +
cval * piplus1[j];
/* need to copy values */
memcpy(piplus1, pi, sizeof(double) * n);
memcpy(pi, pcurr, sizeof(double) * n);
}
/*** now copy the final result ***/
memcpy(result, pcurr, sizeof(double) * n);
}
/* return (pcurr); */
}
BTW, KH: http://www.cs.dartmouth.edu/~geelong/sphere/
|
|
| |
|
|
|
Flip McWho
|
|
New Zealand
May 2002
time: 17:37
|
|
The effort factor is through the roof.
|
|
| |
All times are GMT. The time now is 05:37.
Apolyton Time is 00:37. |
top of page
|
|
|
Forum Rules:
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts
|
HTML code is ON
vB code is ON
Smilies are ON
[IMG] code is ON
|
|
|
|
|
|
08-06-2005 05:53
So it's not an Apolyton problem. Nope. Not at all. No problem here at all. It's you; it's all you, you, you. 
|