23char circhebpimi_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFT991/circhebpimi.C,v 1.4 2014/10/15 12:48:21 j_novak Exp $" ;
146#include "proto_f77.h"
151int* facto_ini(
int ) ;
152double* trigo_ini(
int ) ;
153double* cheb_ini(
const int) ;
154double* chebimp_ini(
const int ) ;
157void circhebpimi(
const int* deg,
const int* dimc,
double* cf,
158 const int* dimf,
double* ff)
177 cout <<
"circhebpimi: nr > n3c : nr = " << nr <<
" , n3c = "
183 cout <<
"circhebpimi: nr > n3f : nr = " << nr <<
" , n3f = "
189 cout <<
"circhebpimi: n1c > n1f : n1c = " << n1c <<
" , n1f = "
195 cout <<
"circhebpimi: n2c > n2f : n2c = " << n2c <<
" , n2f = "
206 int* facto = facto_ini(nm1) ;
207 double* trigo = trigo_ini(nm1) ;
210 double* sinp = cheb_ini(nr);
213 double* x = chebimp_ini(nr);
217 double* g = (
double*)( malloc( (nm1+2)*
sizeof(
double) ) ) ;
218 double* t1 = (
double*)( malloc( (nm1+2)*
sizeof(double) ) ) ;
228 int n2n3f = n2f * n3f ;
229 int n2n3c = n2c * n3c ;
245 for (k=0; k<n2c; k++) {
247 int i0 = n2n3c * j + n3c * k ;
248 double* cf0 = cf + i0 ;
250 i0 = n2n3f * j + n3f * k ;
251 double* ff0 = ff + i0 ;
256 t1[0] = .5 * cf0[0] ;
257 for (i=1; i<nm1; i++) t1[i] = .5 * ( cf0[i] + cf0[i-1] ) ;
258 t1[nm1] = .5 * cf0[nr-2] ;
276 for ( i = 3; i < nr; i += 2 ) {
277 ff0[i] = t1[i] - c1 ;
282 double fmoins0 = nm1s2 * c1 + som ;
288 for ( i = 3; i < nr; i += 2 ) {
289 g[i] = 0.25 * ( ff0[i] - ff0[i-2] ) ;
302 for (i=2; i<nm1; i += 2 ) g[i] = 0.5 * t1[i] ;
309 F77_fft991( g, t1, trigo, facto, &inc, &jump, &nm1, &lot, &isign) ;
314 for ( i = 1; i < nm1s2 ; i++ ) {
320 int ixsym = nm1 - isym ;
322 double fp = .5 * ( g[i] + g[isym] ) ;
323 double fm = .5 * ( g[i] - g[isym] ) / sinp[i] ;
325 ff0[ix] = ( fp + fm ) / x[ix];
326 ff0[ixsym] = ( fp - fm ) / x[ixsym] ;
331 ff0[nm1] = g[0] + fmoins0 ;
332 ff0[nm1s2] = g[nm1s2] / x[nm1s2] ;
342 if ( (j != 1) && (j != n1f-1) ) {
346 for (k=0; k<n2c; k++) {
348 int i0 = n2n3c * j + n3c * k ;
349 double* cf0 = cf + i0 ;
351 i0 = n2n3f * j + n3f * k ;
352 double* ff0 = ff + i0 ;
357 t1[0] = .5 * cf0[0] ;
358 for (i=1; i<nm1; i++) t1[i] = .5 * ( cf0[i] + cf0[i-1] ) ;
359 t1[nm1] = .5 * cf0[nr-2] ;
377 for ( i = 3; i < nr; i += 2 ) {
378 ff0[i] = t1[i] - c1 ;
383 double fmoins0 = nm1s2 * c1 + som ;
389 for ( i = 3; i < nr; i += 2 ) {
390 g[i] = 0.25 * ( ff0[i] - ff0[i-2] ) ;
403 for (i=2; i<nm1; i += 2 ) g[i] = 0.5 * t1[i] ;
410 F77_fft991( g, t1, trigo, facto, &inc, &jump, &nm1, &lot, &isign) ;
415 for ( i = 1; i < nm1s2 ; i++ ) {
421 int ixsym = nm1 - isym ;
423 double fp = .5 * ( g[i] + g[isym] ) ;
424 double fm = .5 * ( g[i] - g[isym] ) / sinp[i] ;
426 ff0[ix] = ( fp + fm ) / x[ix];
427 ff0[ixsym] = ( fp - fm ) / x[ixsym] ;
432 ff0[nm1] = g[0] + fmoins0 ;
433 ff0[nm1s2] = g[nm1s2] / x[nm1s2] ;
463 for (k=0; k<n2c; k++) {
465 int i0 = n2n3c * j + n3c * k ;
466 double* cf0 = cf + i0 ;
468 i0 = n2n3f * j + n3f * k ;
469 double* ff0 = ff + i0 ;
487 for ( i = 3; i < nr; i += 2 ) {
488 ff0[i] = cf0[i] - c1 ;
493 double fmoins0 = nm1s2 * c1 + som ;
499 for ( i = 3; i < nr; i += 2 ) {
500 g[i] = 0.25 * ( ff0[i] - ff0[i-2] ) ;
513 for (i=2; i<nm1; i += 2 ) g[i] = 0.5 * cf0[i] ;
520 F77_fft991( g, t1, trigo, facto, &inc, &jump, &nm1, &lot, &isign) ;
525 for ( i = 1; i < nm1s2 ; i++ ) {
531 int ixsym = nm1 - isym ;
533 double fp = .5 * ( g[i] + g[isym] ) ;
534 double fm = .5 * ( g[i] - g[isym] ) / sinp[i] ;
537 ff0[ixsym] = fp - fm ;
541 ff0[0] = g[0] - fmoins0 ;
542 ff0[nm1] = g[0] + fmoins0 ;
543 ff0[nm1s2] = g[nm1s2] ;
558 for (k=0; k<n2c; k++) {
560 int i0 = n2n3c * j + n3c * k ;
561 double* cf0 = cf + i0 ;
563 i0 = n2n3f * j + n3f * k ;
564 double* ff0 = ff + i0 ;
582 for ( i = 3; i < nr; i += 2 ) {
583 ff0[i] = cf0[i] - c1 ;
588 double fmoins0 = nm1s2 * c1 + som ;
594 for ( i = 3; i < nr; i += 2 ) {
595 g[i] = 0.25 * ( ff0[i] - ff0[i-2] ) ;
608 for (i=2; i<nm1; i += 2 ) g[i] = 0.5 * cf0[i] ;
615 F77_fft991( g, t1, trigo, facto, &inc, &jump, &nm1, &lot, &isign) ;
620 for ( i = 1; i < nm1s2 ; i++ ) {
626 int ixsym = nm1 - isym ;
628 double fp = .5 * ( g[i] + g[isym] ) ;
629 double fm = .5 * ( g[i] - g[isym] ) / sinp[i] ;
632 ff0[ixsym] = fp - fm ;
636 ff0[0] = g[0] - fmoins0 ;
637 ff0[nm1] = g[0] + fmoins0 ;
638 ff0[nm1s2] = g[nm1s2] ;