23char circhebpii_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFTW3/circhebpii.C,v 1.3 2014/10/13 08:53:20 j_novak Exp $" ;
116fftw_plan back_fft(
int, Tbl*&) ;
117double* cheb_ini(
const int) ;
118double* chebimp_ini(
const int ) ;
122void circhebpii(
const int* deg,
const int* dimc,
double* cf,
123 const int* dimf,
double* ff)
142 cout <<
"circhebpii: nr > n3c : nr = " << nr <<
" , n3c = "
148 cout <<
"circhebpii: nr > n3f : nr = " << nr <<
" , n3f = "
154 cout <<
"circhebpii: n1c > n1f : n1c = " << n1c <<
" , n1f = "
160 cout <<
"circhebpii: n2c > n2f : n2c = " << n2c <<
" , n2f = "
172 fftw_plan p = back_fft(nm1, pg) ;
174 double* t1 =
new double[nr] ;
177 double* sinp = cheb_ini(nr);
180 double* x = chebimp_ini(nr);
184 int n2n3f = n2f * n3f ;
185 int n2n3c = n2c * n3c ;
193 int borne_phi = ( n1c > 1 ) ? n1c-1 : 1 ;
195 for (j=0; j< borne_phi; j++) {
202 for (k=1; k<n2c; k+=2) {
204 int i0 = n2n3c * j + n3c * k ;
205 double* cf0 = cf + i0 ;
207 i0 = n2n3f * j + n3f * k ;
208 double* ff0 = ff + i0 ;
222 for ( i = 3; i < nr; i += 2 ) {
223 ff0[i] = cf0[i] - c1 ;
228 double fmoins0 = nm1s2 * c1 + som ;
233 for ( i = 3; i < nr; i += 2 ) {
234 g.set(nm1-i/2) = 0.25 * ( ff0[i] - ff0[i-2] ) ;
246 for (i=1; i<nm1s2; i++) g.set(i) = 0.5 * cf0[2*i] ;
247 g.set(nm1s2) = cf0[nm1] ;
258 for ( i = 1; i < nm1s2 ; i++ ) {
264 int ixsym = nm1 - isym ;
266 double fp = .5 * ( g(i) + g(isym) ) ;
267 double fm = .5 * ( g(i) - g(isym) ) / sinp[i] ;
270 ff0[ixsym] = fp - fm ;
274 ff0[0] = g(0) - fmoins0 ;
275 ff0[nm1] = g(0) + fmoins0 ;
276 ff0[nm1s2] = g(nm1s2) ;
281 for (k=0; k<n2c; k+=2) {
283 int i0 = n2n3c * j + n3c * k ;
284 double* cf0 = cf + i0 ;
286 i0 = n2n3f * j + n3f * k ;
287 double* ff0 = ff + i0 ;
292 t1[0] = .5 * cf0[0] ;
293 for (i=1; i<nm1; i++) t1[i] = .5 * ( cf0[i] + cf0[i-1] ) ;
294 t1[nm1] = .5 * cf0[nr-2] ;
308 for ( i = 3; i < nr; i += 2 ) {
309 ff0[i] = t1[i] - c1 ;
314 double fmoins0 = nm1s2 * c1 + som ;
319 for ( i = 3; i < nr; i += 2 ) {
320 g.set(nm1-i/2) = 0.25 * ( ff0[i] - ff0[i-2] ) ;
332 for (i=1; i<nm1s2; i++ ) g.set(i) = 0.5 * t1[2*i] ;
333 g.set(nm1s2) = t1[nm1] ;
344 for ( i = 1; i < nm1s2 ; i++ ) {
350 int ixsym = nm1 - isym ;
352 double fp = .5 * ( g(i) + g(isym) ) ;
353 double fm = .5 * ( g(i) - g(isym) ) / sinp[i] ;
355 ff0[ix] = ( fp + fm ) / x[ix];
356 ff0[ixsym] = ( fp - fm ) / x[ixsym] ;
361 ff0[nm1] = g(0) + fmoins0 ;
362 ff0[nm1s2] = g(nm1s2) / x[nm1s2] ;