FLAG  1.0b1
Exact Fourier-Laguerre transform in spherical coordinates
flag_spherbessel.c
Go to the documentation of this file.
1 
2 // FLAG package
3 // Copyright (C) 2012
4 // Boris Leistedt & Jason McEwen
5 
6 #include "flag.h"
7 #include <math.h>
8 #include <stdlib.h>
9 #include <complex.h>
10 #include <gsl/gsl_math.h>
11 #include <gsl/gsl_sf.h>
12 
25 void flag_spherlaguerre2spherbessel(double *flk, const double *fn, double *kvalues, int Nk, int N, int ell, double tau)
26 {
27  double PI = 3.141592653589793;
28  double *sbesselslag = (double*)calloc(N*Nk, sizeof(double));
29 
30  flag_sbesselslag(sbesselslag, ell, kvalues, Nk, N, tau);
31 
32  int k, n;
33  for(k = 0; k < Nk; k++){
34  for(n = 0; n < N; n++){
35  flk[k] += sqrt(2.0/PI) * fn[n] * sbesselslag[n * Nk + k];
36  }
37  }
38 
39  free(sbesselslag);
40 
41 }
42 
55 void flag_fourierlaguerre2fourierbessel(complex double *flmk, complex double *flmn, double *kvalues, int Nk, int N, int L, double tau)
56 {
57 
58 }
59 
71 void flag_sbesselslag(double *sbesselslag, int ell, double *kvalues, int Nk, int N, double tau)
72 {
73  int k, n, j, m;
74  long double hold, result, temp, y, fac, corrfac, f_jl2bis;
75  long double cw_jp, c_jp, w_jl, ktilde, a, b, c, z, f_jl2_even, f_jl2_odd, f_jl1_even, f_jl1_odd, f_jl0, f_jl1, f_jl2, T0, T1, T2;
76  long double PI = 3.141592653589793;
77 
78  double *mulk_coefs = (double*)calloc(N+2, sizeof(double));
79  long double *vals = (double*)calloc(N, sizeof(long double));
80 
81  for(k = 0; k < Nk; k++){
82 
83  //flag_mulk(mulk_coefs, N+2, ell, kvalues[k], tau);
84 
85  ktilde = tau * kvalues[k];
86 
87  for(n = 0; n < N; n++){
88 
89 
90  for(j = 0; j <= n; j++){
91 
92  if ( (j % 2) == 0 ){ // j even
93  m = (long double)j / 2.0;
94  a = ((long double)ell + 1) / 2.0;
95  b = (long double)ell / 2.0 + 1;
96  c = (long double)ell + 1.5;
97  z = - 4.0 * pow(ktilde, 2.0);
98  f_jl0 = f_jl1_even;
99  f_jl1 = f_jl2_even;
100  //printf("\n j=%i even ; a=%f, b=%f, c=%f",j,a,b,c);
101  } else { // j odd
102  m = ((long double)j - 1) / 2.0;
103  a = (long double)ell / 2.0 + 1.0;
104  b = ((long double)ell + 3) / 2.0;
105  c = (long double)ell + 1.5;
106  z = - (long double)4 * pow(ktilde, 2.0);
107  f_jl0 = f_jl1_odd;
108  f_jl1 = f_jl2_odd;
109  //printf("\n j=%i odd ; a=%f, b=%f, c=%f",j,a,b,c);
110  }
111 
112  if( j == 0 ){
113  fac = 1.0;
114  c_jp = ((long double)n + 2) * (n + 1) / 2.0 ;
115  w_jl = sqrt(PI) * pow(ktilde, ell) * pow(tau, 3.0) * gsl_sf_fact(ell + 2) ;
116  cw_jp = c_jp * w_jl;
117  corrfac = 1.0;//pow(1-z, -a-1) * gsl_sf_hyperg_2F1_renorm(a+1, c-b-1, c, z/(z-1));
118  f_jl1_even = pow(1-z,-a) * gsl_sf_hyperg_2F1_renorm(c-b,a,c,z/(z-1)) / corrfac;
119  f_jl2_even = pow(1-z, -a-1) * gsl_sf_hyperg_2F1_renorm(a+1, c-b-1, c, z/(z-1)) / corrfac; // pow(1-z,-a-1) * gsl_sf_hyperg_2F1_renorm(c-b-1,a+1,c,z/(z-1));
120  f_jl2 = f_jl2_even;
121  sbesselslag[n * Nk + k] += pow((n+1)*(n+2), -0.5) * corrfac * pow(2,j+2) * fac * cw_jp * f_jl2 ;
122  printf("k=%f n=%i j=%i difference=%Le\n", kvalues[k], n, j, pow((n+1)*(n+2), -0.5) * corrfac * pow(2,j+2) * fac * cw_jp * f_jl2 );
123 
124  }else{
125  fac = ((long double)ell + 2 + j) * (-1) * ((long double)n - j + 1.0) / ((long double)j * (j + 2.0));
126  if( j == 1 ){
127  f_jl2_odd = pow(1-z, -b) * gsl_sf_hyperg_2F1_renorm(b, c-a, c, z/(z-1)) / corrfac; //
128  f_jl1_odd = pow(1-z, -a-1) * gsl_sf_hyperg_2F1_renorm(a+1, c-b-1, c, z/(z-1)) / corrfac;
129  f_jl2bis = f_jl1_odd;
130  sbesselslag[n * Nk + k] += pow((n+1)*(n+2), -0.5) * corrfac * pow(2,j+1) * cw_jp * (f_jl2 + 2 * fac * f_jl2bis);
131  printf("k=%f n=%i j=%i difference=%Le\n", kvalues[k], n, j, pow((n+1)*(n+2), -0.5) * corrfac * pow(2,j+1) * cw_jp * (f_jl2 + 2 * fac * f_jl2bis) );
132  f_jl2 = f_jl2bis;
133  } else {
134  T0 = (c - a - m) * (c - b - m) * (c - a - b - 2*m - 1.0);
135  T1 = (c-a-b-2*m)*( c*(a+b-c+2*m) + c - 2*(a+m)*(b+m) + z*( (a+b+2*m)*(c-a-b-2*m) + 2*(a+m)*(b+m) - c + 1.0 ) );
136  T2 = (a + m) * (b + m) * (c - a - b - 2*m + 1.0) * pow(1 - z, 2.0) ;
137  if ( (j % 2) == 0 ){ // j even
138  f_jl2 = (-T1/T2) * f_jl1 + (-T0/T2) * f_jl0 ;
139  if( j == n ){
140  sbesselslag[n * Nk + k] += pow((n+1)*(n+2), -0.5) * corrfac * pow(2,j+2) * fac * cw_jp * f_jl2 ;
141  printf("k=%f n=%i j=%i difference=%Le\n", kvalues[k], n, j, pow((n+1)*(n+2), -0.5) * corrfac * pow(2,j+2) * fac * cw_jp * f_jl2 );
142  }
143  }else{
144  f_jl2bis = ( (-T1/T2) * f_jl1 + (-T0/T2) * f_jl0 ) ;
145  sbesselslag[n * Nk + k] += pow((n+1)*(n+2), -0.5) * corrfac * pow(2,j+1) * cw_jp * (f_jl2 + 2 * fac * f_jl2bis);
146  printf("k=%f n=%i j=%i difference=%Le\n", kvalues[k], n, j, pow((n+1)*(n+2), -0.5) * pow(2,j+1) * cw_jp * (f_jl2 + 2 * fac * f_jl2bis) );
147  f_jl2 = f_jl2bis;
148  }
149  }
150  cw_jp = fac * cw_jp ;
151  c_jp = - ((long double)n - j + 1.0) / ((long double)j * (j + 2.0)) * c_jp ;
152  w_jl = ((long double)ell + 2 + j) * w_jl;
153  }
154 
155  if ( (j % 2) == 0 ){ // j even
156  f_jl1_even = f_jl2_even;
157  f_jl2_even = f_jl2;
158  } else { // j odd
159  f_jl1_odd = f_jl2_odd;
160  f_jl2_odd = f_jl2;
161  }
162 
163  //vals[j] = pow((n+1)*(n+2), -0.5) * cw_jp * ( pow(2,j+2) * f_jl2 );
164  //printf("\n cw_jp=%e f_jl2=%e val=%e", cw_jp, f_jl2, vals[j]);
165  //printf("\n k=%f n=%i j=%i cw_jp=%Le f_jl2=%Le val=%Le", kvalues[k], n, j, cw_jp, f_jl2, vals[j]);
166  //sbesselslag[n * Nk + k] += vals[j] ;
167 
168  }
169  printf("sbesselslag[n * Nk + k] = %e\n", sbesselslag[n * Nk + k]);
170  }
171  }
172 
173  free(mulk_coefs);
174 
175 }
176 
180 void flag_mulk(double *mulk, int n, int ell, double k, double tau)
181 {
182  double PI = 3.141592653589793;
183  int j;
184  double a, b, fac1, fac2, fac3, rec0, rec1, rec2;
185  double c = (double)ell + 1.5;
186  double ktilde = tau * k;
187  double d = - 4.0 * ktilde * ktilde;
188  double fac = pow(tau, 3.0) * sqrt(PI) * pow(ktilde, ell);
189 
190  // Even
191  a = (ell + 1.0) / 2.0;
192  b = ell / 2.0 + 1.0;
193  for(j = 0; j <= n/2; j++){
194  //printf("j1=%i, 2F1(%f, %f, %f, %f) = ",2*j,a+j,b+j,c,d);fflush(NULL);
195  if( j == 0 ){
196  rec1 = pow(1-d, -b) * gsl_sf_hyperg_2F1_renorm(b, c-a, c, d/(d-1)) ; //* gsl_sf_fact(2*j+ell) * pow(2, 2*j);
197  // rec1 = ( atan(sqrt(-d)) / sqrt(-d) ) / gsl_sf_gamma(c);
198  rec0 = rec1;
199  //printf("%3.3e\n",rec0);fflush(NULL);
200  }else if( j == 1 ){
201  rec0 = pow(1-d, -a-1) * gsl_sf_hyperg_2F1_renorm(a+1, c-b-1, c, d/(d-1)) ; //* gsl_sf_fact(2*j+ell) * pow(2, 2*j);
202  // rec0 = pow(d-1,-2.0) / gsl_sf_gamma(c);
203  //printf("%3.3e\n",rec0);fflush(NULL);
204  }else{
205  rec2 = rec1;
206  rec1 = rec0;
207  fac1 = (c-a-j+1)*(c-b-j+1)*(c-a-b-2*j+1) * rec2 ; //* 16*(2*j+ell)*(2*j+ell-1)*(2*j+ell-2)*(2*j+ell-3);
208  fac2 = (c-a-b-2*j+2)*( c*(a+b-c+2*j-2) + c - 2.0*(a+j-1)*(b+j-1)
209  + d*( (a+b+2*j-2)*(c-a-b-2*j+2) + 2.0*(a+j-1)*(b+j-1) - c + 1 ) ) * rec1 ; //* 4*(2*j+ell)*(2*j+ell-1);
210  fac3 = ( (a+j-1)*(b+j-1)*(c-a-b-2*j+3)*(1.0-d)*(1.0-d) );
211  rec0 = -1.0 * (fac1 + fac2) / fac3 ;
212  //printf("%3.3e\n",rec0);fflush(NULL);
213  //printf("fac1=%5.4f, fac2=%5.4f, fac3=%5.4f\n",fac1,fac2,fac3);fflush(NULL);
214  }
215  mulk[2*j] = rec0; //fac * //* gsl_sf_gamma(2*j+ell+1) * pow(2, 2*j)
216  }
217 
218  // Odd
219  a = ell / 2.0 + 1.0;
220  b = (ell + 3) / 2.0;
221  for(j = 0; j <= n/2-1; j++){
222  //printf("j2=%i, 2F1(%f, %f, %f, %f) = ",2*j+1,a+j,b+j,c,d);fflush(NULL);
223  if( j == 0 ){
224  rec1 = pow(1-d, -b) * gsl_sf_hyperg_2F1_renorm(b, c-a, c, d/(d-1)) ; //* gsl_sf_fact(2*j+ell+1) * pow(2, 2*j+1);
225  //rec1 = pow(1-d, -1.0) / gsl_sf_gamma(c) ;
226  rec0 = rec1;
227  //printf("%3.3e\n",rec0);fflush(NULL);
228  }else if( j == 1 ){
229  rec0 = pow(1-d, -a-1) * gsl_sf_hyperg_2F1_renorm(a+1, c-b-1, c, d/(d-1)) ; //* gsl_sf_fact(2*j+ell+1) * pow(2, 2*j+1);
230  //rec0 = ((d+3)/(3*pow(1-d, 3.0))) / gsl_sf_gamma(c);//
231  //printf("%3.3e\n",rec0);fflush(NULL);
232  }else{
233  rec2 = rec1;
234  rec1 = rec0;
235  fac1 = (c-a-j+1)*(c-b-j+1)*(c-a-b-2*j+1) * rec2 ; //* 16*(2*j+ell)*(2*j+ell-1)*(2*j+ell-2)*(2*j+ell-3);
236  fac2 = (c-a-b-2*j+2)*( c*(a+b-c+2*j-2) + c - 2.0*(a+j-1)*(b+j-1)
237  + d*( (a+b+2*j-2)*(c-a-b-2*j+2) + 2.0*(a+j-1)*(b+j-1) - c + 1 ) ) * rec1 ; //* 4*(2*j+ell)*(2*j+ell-1);
238  fac3 = ( (a+j-1)*(b+j-1)*(c-a-b-2*j+3)*(1.0-d)*(1.0-d) );
239  rec0 = -1.0 * (fac1 + fac2) / fac3 ;
240  //printf("%3.3e\n",rec0);fflush(NULL);
241  //printf("fac1=%5.4f, fac2=%5.4f, fac3=%5.4f\n",fac1,fac2,fac3);fflush(NULL);
242  }
243  mulk[2*j+1] = rec0;//fac * * gsl_sf_gamma(2*j+ell+2) * pow(2, 2*j+1)
244  }
245 
246  for(j = 0; j <= n; j++){
247  //printf("j=%i, mulk(%1.1f, %1.1f, %1.1f, %1.1f) = %6.5e\n",j,(ell+j+1)/2.0,(ell+j)/2.0+1,c,d,mulk[j]);fflush(NULL);
248  }
249 
250 }
251 
252 
253 void bubbleSort(double numbers[], int array_size)
254 {
255  int i, j;
256  double temp;
257 
258  for (i = (array_size - 1); i > 0; i--)
259  {
260  for (j = 1; j <= i; j++)
261  {
262  if ( numbers[j-1] > numbers[j] )
263  {
264  temp = numbers[j-1];
265  numbers[j-1] = numbers[j];
266  numbers[j] = temp;
267  }
268  }
269  }
270 }
271 
272 
273 void flag_sbesselslag_backup(double *sbesselslag, int ell, double *kvalues, int Nk, int N, double tau)
274 {
275  int k, n, j;
276  double weight, hold, result, temp, y;
277 
278  double *mulk_coefs = (double*)calloc(N+2, sizeof(double));
279  double *vals = (double*)calloc(N, sizeof(double));
280 
281  for(k = 0; k < Nk; k++){
282  flag_mulk(mulk_coefs, N+2, ell, kvalues[k], tau);
283  for(n = 0; n < N; n++){
284  printf("\n");
285  for(j = 0; j <= n; j++){
286  if( j == 0 ){
287  weight = pow(2, 2) * gsl_sf_fact(ell + 2) * (n + 1.0) * (n + 2.0) / 2.0 ; //
288  }else{
289  weight = - 2 * (ell + 2 + j) * (n - j + 1) * weight / (j * (j + 2)); //
290  }
291  sbesselslag[n * Nk + k] += pow((n+1)*(n+2), -0.5) * weight * mulk_coefs[j+2] ;
292  printf(" %3.0e ",weight);
293  vals[j] = pow((n+1)*(n+2), -0.5) * ( weight * mulk_coefs[j+2] );
294  }
295  result = 0.0, hold = 0.0, temp = 0.0, y = 0.0;
296  for(j = 0; j <= n/2; j++){
297  if( 2*j+1 <= n ){
298  y = (vals[2*j] + vals[2*j+1]) - hold;
299  temp = result + y;
300  hold = (temp - result) - y;
301  result = temp;
302  } else {
303  y = vals[2*j] - hold;
304  temp = result + y;
305  hold = (temp - result) - y;
306  result = temp;
307  }
308  // printf("n = %i, j = %i, val = %3.2e, val = %3.2e, result = ", n, j, vals[2*j], vals[2*j+1]);
309  // printf("%3.2e\n",result);
310  }
311  /*
312  bubbleSort(vals, n+1);
313  for(j = 0; j <= n; j++)
314  printf(" %3.2e \n", vals[j]);
315  result = 0.0, hold = 0.0, temp = 0.0, y = 0.0;
316  for(j = 0; j <= n/2; j++){
317  y = (vals[j] + vals[n-j]) - hold;
318  temp = result + y;
319  hold = (temp - result) - y;
320  result = temp;
321  printf("n = %i, j = %i, val = %3.2e, val = %3.2e, result = ", n, j, vals[j], vals[n-j]);
322  printf("%3.2e\n",result);
323  }
324  if( n % 2 != 1 ){
325  y = - vals[n/2] - hold;
326  temp = result + y;
327  hold = (temp - result) - y;
328  result = temp;
329  }
330  */
331  sbesselslag[n * Nk + k] = result ;
332  }
333  }
334 
335  free(mulk_coefs);
336 
337 }
338 
339 
351 void flag_sbesselslag_backup2(double *sbesselslag, int ell, double *kvalues, int Nk, int N, double tau)
352 {
353  int k, n, j, m;
354  double hold, result, temp, y;
355  double c_jp, w_jl, ktilde, a, b, c, z, f_jl2_even, f_jl2_odd, f_jl1_even, f_jl1_odd, f_jl0, f_jl1, f_jl2, T0, T1, T2;
356  double PI = 3.141592653589793;
357 
358  double *mulk_coefs = (double*)calloc(N+2, sizeof(double));
359  double *vals = (double*)calloc(N, sizeof(double));
360 
361  for(k = 0; k < Nk; k++){
362 
363  flag_mulk(mulk_coefs, N+2, ell, kvalues[k], tau);
364 
365  ktilde = tau * kvalues[k];
366 
367  for(n = 0; n < N; n++){
368 
369  printf("\n n = %i", n);
370 
371  for(j = 0; j <= n; j++){
372 
373  if ( (j % 2) == 0 ){ // j even
374  m = j / 2.0;
375  a = (ell + 1) / 2.0;
376  b = ell / 2.0 + 1;
377  c = ell + 1.5;
378  z = - 4.0 * pow(ktilde, 2.0);
379  f_jl0 = f_jl1_even;
380  f_jl1 = f_jl2_even;
381  //printf("\n j=%i even ; a=%f, b=%f, c=%f",j,a,b,c);
382  } else { // j odd
383  m = ((double)j - 1) / 2.0;
384  a = (double)ell / 2.0 + 1.0;
385  b = ((double)ell + 3) / 2.0;
386  c = (double)ell + 1.5;
387  z = - (double)4 * pow(ktilde, 2.0);
388  f_jl0 = f_jl1_odd;
389  f_jl1 = f_jl2_odd;
390  //printf("\n j=%i odd ; a=%f, b=%f, c=%f",j,a,b,c);
391  }
392 
393  if( j == 0 ){
394  c_jp = (n + 2) * (n + 1) / 2.0 ;
395  w_jl = 4 * sqrt(PI) * pow(ktilde, ell) * pow(tau, 1.5) * gsl_sf_fact(ell + 2) ;
396  f_jl1_even = pow(1-z, -b) * gsl_sf_hyperg_2F1_renorm(b, c-a, c, z/(z-1)); // pow(1-z,-a) * gsl_sf_hyperg_2F1_renorm(c-b,a,c,z/(z-1));
397  f_jl2_even = pow(1-z, -a-1) * gsl_sf_hyperg_2F1_renorm(a+1, c-b-1, c, z/(z-1)); // pow(1-z,-a-1) * gsl_sf_hyperg_2F1_renorm(c-b-1,a+1,c,z/(z-1));
398  f_jl2 = f_jl2_even;
399  }else{
400  c_jp = - (n - j + 1) * c_jp / (j * (j + 2));
401  w_jl = 2 * (ell + 2 + j) * w_jl;
402  if( j == 1 ){
403  f_jl2_odd = pow(1-z, -b) * gsl_sf_hyperg_2F1_renorm(b, c-a, c, z/(z-1)); // pow(1-z,-a) * gsl_sf_hyperg_2F1_renorm(c-b,a,c,z/(z-1));
404  f_jl1_odd = pow(1-z, -a-1) * gsl_sf_hyperg_2F1_renorm(a+1, c-b-1, c, z/(z-1)) ; //* gsl_sf_fact(2*j+ell+1) * pow(2, 2*j+1);
405  //f_jl2_odd = pow(1-z, -a-1) * gsl_sf_hyperg_2F1_renorm(a+1, c-b-1, c, z/(z-1)); // pow(1-z,-a-1) * gsl_sf_hyperg_2F1_renorm(c-b-1,a+1,c,z/(z-1));
406  f_jl2 = f_jl1_odd;
407  } else {
408  T0 = (c - a - m) * (c - b - m) * (c - a - b - 2*m - 1.0);
409  T1 = (c-a-b-2*m)*( c*(a+b-c+2*m) + c - 2*(a+m)*(b+m) + z*( (a+b+2*m)*(c-a-b-2*m) + 2*(a+m)*(b+m) - c + 1.0 ) );
410  T2 = (a + m) * (b + m) * (c - a - b - 2*m + 1.0) * pow(1 - z, 2.0) ;
411  f_jl2 = (-T1/T2) * f_jl1 + (-T0/T2) * f_jl0 ;
412  }
413 
414  if ( (j % 2) == 0 ){ // j even
415  f_jl1_even = f_jl2_even;
416  f_jl2_even = f_jl2;
417  } else { // j odd
418  f_jl1_odd = f_jl2_odd;
419  f_jl2_odd = f_jl2;
420  }
421 
422  vals[j] = c_jp * w_jl * f_jl2 ;
423  printf("\n c_jp = %e w_jl = %e f_jl2 = %e mu = %e val = %f", c_jp, w_jl, f_jl2);
424  sbesselslag[n * Nk + k] += vals[j];
425  sbesselslag[n * Nk + k] *= pow((n+1)*(n+2), -0.5);
426  }
427 
428  }
429  /*
430  result = 0.0, hold = 0.0, temp = 0.0, y = 0.0;
431  for(j = 0; j <= n/2; j++){
432  if( 2*j+1 <= n ){
433  y = (vals[2*j] + vals[2*j+1]) - hold;
434  temp = result + y;
435  hold = (temp - result) - y;
436  result = temp;
437  } else {
438  y = vals[2*j] - hold;
439  temp = result + y;
440  hold = (temp - result) - y;
441  result = temp;
442  }
443  // printf("n = %i, j = %i, val = %3.2e, val = %3.2e, result = ", n, j, vals[2*j], vals[2*j+1]);
444  // printf("%3.2e\n",result);
445  }
446  */
447  //sbesselslag[n * Nk + k] = result ;
448  }
449  }
450 
451  free(mulk_coefs);
452 
453 }
454 
455 
456 
457 void flag_sbesselslag_backup3(double *sbesselslag, int ell, double *kvalues, int Nk, int N, double tau)
458 {
459  int k, n, j, m;
460  long double hold, result, temp, y, fac, corrfac, f_jl2bis;
461  long double cw_jp, c_jp, w_jl, ktilde, a, b, c, z, f_jl2_even, f_jl2_odd, f_jl1_even, f_jl1_odd, f_jl0, f_jl1, f_jl2, T0, T1, T2;
462  long double PI = 3.141592653589793;
463 
464  double *mulk_coefs = (double*)calloc(N+2, sizeof(double));
465  long double *vals = (double*)calloc(N, sizeof(long double));
466 
467  for(k = 0; k < Nk; k++){
468 
469  //flag_mulk(mulk_coefs, N+2, ell, kvalues[k], tau);
470 
471  ktilde = tau * kvalues[k];
472 
473  for(n = 0; n < N; n++){
474 
475  //printf("\n n = %i", n);
476 
477  for(j = 0; j <= n; j++){
478 
479  if ( (j % 2) == 0 ){ // j even
480  m = (long double)j / 2.0;
481  a = ((long double)ell + 1) / 2.0;
482  b = (long double)ell / 2.0 + 1;
483  c = (long double)ell + 1.5;
484  z = - 4.0 * pow(ktilde, 2.0);
485  f_jl0 = f_jl1_even;
486  f_jl1 = f_jl2_even;
487  //printf("\n j=%i even ; a=%f, b=%f, c=%f",j,a,b,c);
488  } else { // j odd
489  m = ((long double)j - 1) / 2.0;
490  a = (long double)ell / 2.0 + 1.0;
491  b = ((long double)ell + 3) / 2.0;
492  c = (long double)ell + 1.5;
493  z = - (long double)4 * pow(ktilde, 2.0);
494  f_jl0 = f_jl1_odd;
495  f_jl1 = f_jl2_odd;
496  //printf("\n j=%i odd ; a=%f, b=%f, c=%f",j,a,b,c);
497  }
498 
499  if( j == 0 ){
500  fac = 1.0;
501  c_jp = ((long double)n + 2) * (n + 1) / 2.0 ;
502  w_jl = sqrt(PI) * pow(ktilde, ell) * pow(tau, 3.0) * gsl_sf_fact(ell + 2) ;
503  cw_jp = c_jp * w_jl;
504  corrfac = 1.0;//pow(1-z, -a-1) * gsl_sf_hyperg_2F1_renorm(a+1, c-b-1, c, z/(z-1));
505  f_jl1_even = pow(1-z,-a) * gsl_sf_hyperg_2F1_renorm(c-b,a,c,z/(z-1)) / corrfac;
506  f_jl2_even = pow(1-z, -a-1) * gsl_sf_hyperg_2F1_renorm(a+1, c-b-1, c, z/(z-1)) / corrfac; // pow(1-z,-a-1) * gsl_sf_hyperg_2F1_renorm(c-b-1,a+1,c,z/(z-1));
507  f_jl2 = f_jl2_even;
508  }else{
509  fac = ((long double)ell + 2 + j) * (-1) * ((long double)n - j + 1.0) / ((long double)j * (j + 2.0));
510  if( j == 1 ){
511  f_jl2_odd = pow(1-z, -b) * gsl_sf_hyperg_2F1_renorm(b, c-a, c, z/(z-1)) / corrfac; //
512  f_jl1_odd = pow(1-z, -a-1) * gsl_sf_hyperg_2F1_renorm(a+1, c-b-1, c, z/(z-1)) / corrfac;
513  f_jl2bis = f_jl1_odd;
514  sbesselslag[n * Nk + k] += pow((n+1)*(n+2), -0.5) * corrfac * pow(2,j+1) * cw_jp * (f_jl2 + 2 * fac * f_jl2bis);
515  f_jl2 = f_jl2bis;
516  } else {
517  T0 = (c - a - m) * (c - b - m) * (c - a - b - 2*m - 1.0);
518  T1 = (c-a-b-2*m)*( c*(a+b-c+2*m) + c - 2*(a+m)*(b+m) + z*( (a+b+2*m)*(c-a-b-2*m) + 2*(a+m)*(b+m) - c + 1.0 ) );
519  T2 = (a + m) * (b + m) * (c - a - b - 2*m + 1.0) * pow(1 - z, 2.0) ;
520  if ( (j % 2) == 0 ){ // j even
521  f_jl2 = (-T1/T2) * f_jl1 + (-T0/T2) * f_jl0 ;
522  if( j == n ){
523  sbesselslag[n * Nk + k] += pow((n+1)*(n+2), -0.5) * corrfac * pow(2,j+2) * fac * cw_jp * f_jl2 ;
524  }
525  }else{
526  f_jl2bis = ( (-T1/T2) * f_jl1 + (-T0/T2) * f_jl0 ) ;
527  sbesselslag[n * Nk + k] += pow((n+1)*(n+2), -0.5) * corrfac * pow(2,j+1) * cw_jp * (f_jl2 + 2 * fac * f_jl2bis);
528  printf("n=%i j=%i difference=%Le\n", n, j, pow((n+1)*(n+2), -0.5) * pow(2,j+1) * cw_jp * (f_jl2 + 2 * fac * f_jl2bis) );
529  f_jl2 = f_jl2bis;
530  }
531  }
532  cw_jp = fac * cw_jp ;
533  c_jp = - ((long double)n - j + 1.0) / ((long double)j * (j + 2.0)) * c_jp ;
534  w_jl = ((long double)ell + 2 + j) * w_jl;
535  }
536 
537  if ( (j % 2) == 0 ){ // j even
538  f_jl1_even = f_jl2_even;
539  f_jl2_even = f_jl2;
540  } else { // j odd
541  f_jl1_odd = f_jl2_odd;
542  f_jl2_odd = f_jl2;
543  }
544 
545  //vals[j] = pow((n+1)*(n+2), -0.5) * cw_jp * ( pow(2,j+2) * f_jl2 );
546  //printf("\n cw_jp=%e f_jl2=%e val=%e", cw_jp, f_jl2, vals[j]);
547  //printf("\n k=%f n=%i j=%i cw_jp=%Le f_jl2=%Le val=%Le", kvalues[k], n, j, cw_jp, f_jl2, vals[j]);
548  //sbesselslag[n * Nk + k] += vals[j] ;
549 
550 
551  }
552  printf("sbesselslag[n * Nk + k] = %e\n", sbesselslag[n * Nk + k]);
553  }
554  }
555 
556  free(mulk_coefs);
557 
558 }
flag_sbesselslag
void flag_sbesselslag(double *sbesselslag, int ell, double *kvalues, int Nk, int N, double tau)
Definition: flag_spherbessel.c:71
flag_mulk
void flag_mulk(double *mulk, int n, int ell, double k, double tau)
Definition: flag_spherbessel.c:180
flag_sbesselslag_backup2
void flag_sbesselslag_backup2(double *sbesselslag, int ell, double *kvalues, int Nk, int N, double tau)
Definition: flag_spherbessel.c:351
flag_spherlaguerre2spherbessel
void flag_spherlaguerre2spherbessel(double *flk, const double *fn, double *kvalues, int Nk, int N, int ell, double tau)
Definition: flag_spherbessel.c:25
flag_fourierlaguerre2fourierbessel
void flag_fourierlaguerre2fourierbessel(complex double *flmk, complex double *flmn, double *kvalues, int Nk, int N, int L, double tau)
Definition: flag_spherbessel.c:55
flag_sbesselslag_backup3
void flag_sbesselslag_backup3(double *sbesselslag, int ell, double *kvalues, int Nk, int N, double tau)
Definition: flag_spherbessel.c:457
flag_sbesselslag_backup
void flag_sbesselslag_backup(double *sbesselslag, int ell, double *kvalues, int Nk, int N, double tau)
Definition: flag_spherbessel.c:273
bubbleSort
void bubbleSort(double numbers[], int array_size)
Definition: flag_spherbessel.c:253
Generated on Tue Sep 22 2015 17:01:04 for FLAG by   doxygen 1.8.10