23char circhebpip_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFT991/circhebpip.C,v 1.3 2014/10/13 08:53:17 j_novak Exp $" ;
114#include "proto_f77.h"
118int* facto_ini(
int ) ;
119double* trigo_ini(
int ) ;
120double* cheb_ini(
const int) ;
121double* chebimp_ini(
const int ) ;
124void circhebpip(
const int* deg,
const int* dimc,
double* cf,
125 const int* dimf,
double* ff)
144 cout <<
"circhebpip: nr > n3c : nr = " << nr <<
" , n3c = "
150 cout <<
"circhebpip: nr > n3f : nr = " << nr <<
" , n3f = "
156 cout <<
"circhebpip: n1c > n1f : n1c = " << n1c <<
" , n1f = "
162 cout <<
"circhebpip: n2c > n2f : n2c = " << n2c <<
" , n2f = "
173 int* facto = facto_ini(nm1) ;
174 double* trigo = trigo_ini(nm1) ;
177 double* sinp = cheb_ini(nr);
180 double* x = chebimp_ini(nr);
184 double* g = (
double*)( malloc( (nm1+2)*
sizeof(
double) ) ) ;
185 double* t1 = (
double*)( malloc( (nm1+2)*
sizeof(double) ) ) ;
195 int n2n3f = n2f * n3f ;
196 int n2n3c = n2c * n3c ;
204 int borne_phi = ( n1c > 1 ) ? n1c-1 : 1 ;
206 for (j=0; j< borne_phi; j++) {
213 for (k=0; k<n2c; k+=2) {
215 int i0 = n2n3c * j + n3c * k ;
216 double* cf0 = cf + i0 ;
218 i0 = n2n3f * j + n3f * k ;
219 double* ff0 = ff + i0 ;
233 for ( i = 3; i < nr; i += 2 ) {
234 ff0[i] = cf0[i] - c1 ;
239 double fmoins0 = nm1s2 * c1 + som ;
245 for ( i = 3; i < nr; i += 2 ) {
246 g[i] = 0.25 * ( ff0[i] - ff0[i-2] ) ;
259 for (i=2; i<nm1; i += 2 ) g[i] = 0.5 * cf0[i] ;
266 F77_fft991( g, t1, trigo, facto, &inc, &jump, &nm1, &lot, &isign) ;
271 for ( i = 1; i < nm1s2 ; i++ ) {
277 int ixsym = nm1 - isym ;
279 double fp = .5 * ( g[i] + g[isym] ) ;
280 double fm = .5 * ( g[i] - g[isym] ) / sinp[i] ;
283 ff0[ixsym] = fp - fm ;
287 ff0[0] = g[0] - fmoins0 ;
288 ff0[nm1] = g[0] + fmoins0 ;
289 ff0[nm1s2] = g[nm1s2] ;
295 for (k=1; k<n2c; k+=2) {
297 int i0 = n2n3c * j + n3c * k ;
298 double* cf0 = cf + i0 ;
300 i0 = n2n3f * j + n3f * k ;
301 double* ff0 = ff + i0 ;
306 t1[0] = .5 * cf0[0] ;
307 for (i=1; i<nm1; i++) t1[i] = .5 * ( cf0[i] + cf0[i-1] ) ;
308 t1[nm1] = .5 * cf0[nr-2] ;
322 for ( i = 3; i < nr; i += 2 ) {
323 ff0[i] = t1[i] - c1 ;
328 double fmoins0 = nm1s2 * c1 + som ;
334 for ( i = 3; i < nr; i += 2 ) {
335 g[i] = 0.25 * ( ff0[i] - ff0[i-2] ) ;
348 for (i=2; i<nm1; i += 2 ) g[i] = 0.5 * t1[i] ;
355 F77_fft991( g, t1, trigo, facto, &inc, &jump, &nm1, &lot, &isign) ;
360 for ( i = 1; i < nm1s2 ; i++ ) {
366 int ixsym = nm1 - isym ;
368 double fp = .5 * ( g[i] + g[isym] ) ;
369 double fm = .5 * ( g[i] - g[isym] ) / sinp[i] ;
371 ff0[ix] = ( fp + fm ) / x[ix];
372 ff0[ixsym] = ( fp - fm ) / x[ixsym] ;
377 ff0[nm1] = g[0] + fmoins0 ;
378 ff0[nm1s2] = g[nm1s2] / x[nm1s2] ;