23char circhebpimp_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFTW3/circhebpimp.C,v 1.4 2014/10/13 08:53:20 j_novak Exp $" ;
149fftw_plan back_fft(
int, Tbl*&) ;
150double* cheb_ini(
const int) ;
151double* chebimp_ini(
const int ) ;
154void circhebpimp(
const int* deg,
const int* dimc,
double* cf,
155 const int* dimf,
double* ff)
173 cout <<
"circhebpimp: nr > n3c : nr = " << nr <<
" , n3c = "
179 cout <<
"circhebpimp: nr > n3f : nr = " << nr <<
" , n3f = "
185 cout <<
"circhebpimp: n1c > n1f : n1c = " << n1c <<
" , n1f = "
191 cout <<
"circhebpimp: n2c > n2f : n2c = " << n2c <<
" , n2f = "
203 fftw_plan p = back_fft(nm1, pg) ;
205 double* t1 =
new double[nr] ;
208 double* sinp = cheb_ini(nr);
211 double* x = chebimp_ini(nr);
215 int n2n3f = n2f * n3f ;
216 int n2n3c = n2c * n3c ;
231 for (k=0; k<n2c; k++) {
233 int i0 = n2n3c * j + n3c * k ;
234 double* cf0 = cf + i0 ;
236 i0 = n2n3f * j + n3f * k ;
237 double* ff0 = ff + i0 ;
255 for ( i = 3; i < nr; i += 2 ) {
256 ff0[i] = cf0[i] - c1 ;
261 double fmoins0 = nm1s2 * c1 + som ;
266 for ( i = 3; i < nr; i += 2 ) {
267 g.set(nm1-i/2) = 0.25 * ( ff0[i] - ff0[i-2] ) ;
279 for (i=1; i<nm1s2; i++) g.set(i) = 0.5 * cf0[2*i] ;
280 g.set(nm1s2) = cf0[nm1] ;
291 for ( i = 1; i < nm1s2 ; i++ ) {
297 int ixsym = nm1 - isym ;
299 double fp = .5 * ( g(i) + g(isym) ) ;
300 double fm = .5 * ( g(i) - g(isym) ) / sinp[i] ;
303 ff0[ixsym] = fp - fm ;
307 ff0[0] = g(0) - fmoins0 ;
308 ff0[nm1] = g(0) + fmoins0 ;
309 ff0[nm1s2] = g(nm1s2) ;
319 if ( (j != 1) && (j != n1f-1) ) {
323 for (k=0; k<n2c; k++) {
325 int i0 = n2n3c * j + n3c * k ;
326 double* cf0 = cf + i0 ;
328 i0 = n2n3f * j + n3f * k ;
329 double* ff0 = ff + i0 ;
347 for ( i = 3; i < nr; i += 2 ) {
348 ff0[i] = cf0[i] - c1 ;
353 double fmoins0 = nm1s2 * c1 + som ;
358 for ( i = 3; i < nr; i += 2 ) {
359 g.set(nm1-i/2) = 0.25 * ( ff0[i] - ff0[i-2] ) ;
371 for (i=1; i<nm1s2; i++) g.set(i) = 0.5 * cf0[2*i] ;
372 g.set(nm1s2) = cf0[nm1] ;
383 for ( i = 1; i < nm1s2 ; i++ ) {
389 int ixsym = nm1 - isym ;
391 double fp = .5 * ( g(i) + g(isym) ) ;
392 double fm = .5 * ( g(i) - g(isym) ) / sinp[i] ;
395 ff0[ixsym] = fp - fm ;
399 ff0[0] = g(0) - fmoins0 ;
400 ff0[nm1] = g(0) + fmoins0 ;
401 ff0[nm1s2] = g(nm1s2) ;
431 for (k=0; k<n2c; k++) {
433 int i0 = n2n3c * j + n3c * k ;
434 double* cf0 = cf + i0 ;
436 i0 = n2n3f * j + n3f * k ;
437 double* ff0 = ff + i0 ;
443 t1[0] = .5 * cf0[0] ;
444 for (i=1; i<nm1; i++) t1[i] = .5 * ( cf0[i] + cf0[i-1] ) ;
445 t1[nm1] = .5 * cf0[nr-2] ;
463 for ( i = 3; i < nr; i += 2 ) {
464 ff0[i] = t1[i] - c1 ;
469 double fmoins0 = nm1s2 * c1 + som ;
474 for ( i = 3; i < nr; i += 2 ) {
475 g.set(nm1-i/2) = 0.25 * ( ff0[i] - ff0[i-2] ) ;
487 for (i=1; i<nm1s2; i++ ) g.set(i) = 0.5 * t1[2*i] ;
488 g.set(nm1s2) = t1[nm1] ;
499 for ( i = 1; i < nm1s2 ; i++ ) {
505 int ixsym = nm1 - isym ;
507 double fp = .5 * ( g(i) + g(isym) ) ;
508 double fm = .5 * ( g(i) - g(isym) ) / sinp[i] ;
510 ff0[ix] = ( fp + fm ) / x[ix];
511 ff0[ixsym] = ( fp - fm ) / x[ixsym] ;
516 ff0[nm1] = g(0) + fmoins0 ;
517 ff0[nm1s2] = g(nm1s2) / x[nm1s2] ;
530 for (k=0; k<n2c; k++) {
532 int i0 = n2n3c * j + n3c * k ;
533 double* cf0 = cf + i0 ;
535 i0 = n2n3f * j + n3f * k ;
536 double* ff0 = ff + i0 ;
541 t1[0] = .5 * cf0[0] ;
542 for (i=1; i<nm1; i++) t1[i] = .5 * ( cf0[i] + cf0[i-1] ) ;
543 t1[nm1] = .5 * cf0[nr-2] ;
561 for ( i = 3; i < nr; i += 2 ) {
562 ff0[i] = t1[i] - c1 ;
567 double fmoins0 = nm1s2 * c1 + som ;
572 for ( i = 3; i < nr; i += 2 ) {
573 g.set(nm1-i/2) = 0.25 * ( ff0[i] - ff0[i-2] ) ;
585 for (i=1; i<nm1s2; i++ ) g.set(i) = 0.5 * t1[2*i] ;
586 g.set(nm1s2) = t1[nm1] ;
597 for ( i = 1; i < nm1s2 ; i++ ) {
603 int ixsym = nm1 - isym ;
605 double fp = .5 * ( g(i) + g(isym) ) ;
606 double fm = .5 * ( g(i) - g(isym) ) / sinp[i] ;
608 ff0[ix] = ( fp + fm ) / x[ix];
609 ff0[ixsym] = ( fp - fm ) / x[ixsym] ;
614 ff0[nm1] = g(0) + fmoins0 ;
615 ff0[nm1s2] = g(nm1s2) / x[nm1s2] ;