FLAG  1.0b1
Exact Fourier-Laguerre transform in spherical coordinates
flag_core.c
Go to the documentation of this file.
1 // FLAG package
2 // Copyright (C) 2012
3 // Boris Leistedt & Jason McEwen
4 
5 #include "flag.h"
6 #include <math.h>
7 #include <stdlib.h>
8 #include <complex.h>
9 #include <fftw3.h>
10 #include <ssht.h>
11 #include <assert.h>
12 
13 int ssht_fr_size_mw(int L)
14 {// In case we want to extend to various sampling schemes.
15  assert(L > 0);
16  int mapsize = L*(2*L-1); // MW sampling scheme
17  return mapsize;
18 }
19 
20 int ssht_flm_size(int L)
21 {// In case we want to extend to various sampling schemes
22  assert(L > 0);
23  int mapsize = L*L;
24  return mapsize;
25 }
26 
27 int flag_core_flmn_size(int L, int N)
28 {// In case we want to extend to various sampling schemes
29  assert(L > 0);
30  assert(N > 1);
31  return ssht_flm_size(L)*N;
32 }
33 
34 int flag_core_f_size_mw(int L, int N)
35 {// In case we want to extend to various sampling schemes
36  assert(L > 0);
37  assert(N > 1);
38  return ssht_fr_size_mw(L) * (N);
39 }
40 
41 void flag_core_allocate_flmn(complex double **flmn, int L, int N)
42 {
43  assert(L > 0);
44  assert(N > 1);
45  int flmsize = ssht_flm_size(L);
46  long totalsize = N*flmsize;
47  *flmn = (complex double*)calloc(totalsize, sizeof(complex double));
48  assert(flmn != NULL);
49 }
50 
51 void flag_core_allocate_f(complex double **f, int L, int N)
52 {
53  assert(L > 0);
54  assert(N > 1);
55  int frsize = ssht_fr_size_mw(L);
56  long totalsize = (N) * frsize;
57  *f = (complex double*)calloc(totalsize, sizeof(complex double));
58  assert(f != NULL);
59 }
60 
61 void flag_core_allocate_f_real(double **f, int L, int N)
62 {
63  assert(L > 0);
64  assert(N > 1);
65  int frsize = ssht_fr_size_mw(L);
66  long totalsize = (N) * frsize;
67  *f = (double*)calloc(totalsize, sizeof(double));
68  assert(f != NULL);
69 }
70 
71 void flag_core_analysis(complex double *flmn,
72  const complex double *f,
73  int L, double tau, int N, int spin)
74 {
75  assert(L > 0);
76  assert(N > 1);
77  //const int alpha = ALPHA;
78  int verbosity = 0;
79  int n;
80  int flmsize = ssht_flm_size(L);
81  int frsize = ssht_fr_size_mw(L);
82  ssht_dl_method_t dl_method = SSHT_DL_TRAPANI;
83  int offset_lm, offset_r;
84 
85  complex double *flmr;
86  flag_core_allocate_flmn(&flmr, L, N);
87 
88  for (n = 0; n < N; n++){
89  //printf("> Analysis: layer %i on %i\n", n+1,N+1);
90  offset_lm = n * flmsize;
91  offset_r = n * frsize;
92  ssht_core_mw_forward_sov_conv_sym(flmr + offset_lm, f + offset_r, L, spin, dl_method, verbosity);
93  }
94 
95  double *nodes = (double*)calloc(N, sizeof(double));
96  double *weights = (double*)calloc(N, sizeof(double));
97  assert(nodes != NULL);
98  assert(weights != NULL);
99  flag_spherlaguerre_sampling(nodes, weights, tau, N);
100  //printf("> Mapped spherical Laguerre transform...");
101  fflush(NULL);
102  flag_spherlaguerre_mapped_analysis(flmn, flmr, nodes, weights, tau, N, flmsize);
103  //printf("done\n");
104  free(nodes);
105  free(weights);
106  free(flmr);
107 
108 }
109 
110 void flag_core_synthesis(complex double *f,
111  const complex double *flmn,
112  const double *nodes, int Nnodes,
113  int L, double tau, int N, int spin)
114 {
115  assert(L > 0);
116  assert(N > 1);
117  assert(nodes != NULL);
118  //const int alpha = ALPHA;
119  int verbosity = 0;
120  int n, offset_lm, offset_r;
121  int flmsize = ssht_flm_size(L);
122  int frsize = ssht_fr_size_mw(L);
123  ssht_dl_method_t dl_method = SSHT_DL_TRAPANI;
124 
125  complex double *flmr;
126  flag_core_allocate_flmn(&flmr, L, Nnodes);
127  //printf("> Mapped spherical Laguerre transform...");fflush(NULL);
128  flag_spherlaguerre_mapped_synthesis(flmr, flmn, nodes, Nnodes, tau, N, flmsize);
129  //printf("done\n");
130 
131  for (n = 0; n < Nnodes; n++){
132  //printf("> Synthesis: layer %i on %i\n",n+1,Nnodes);
133  offset_lm = n * flmsize;
134  offset_r = n * frsize;
135  ssht_core_mw_inverse_sov_sym(f + offset_r, flmr + offset_lm, L, spin, dl_method, verbosity);
136  }
137 
138  free(flmr);
139 }
140 
141 void flag_core_analysis_real(complex double *flmn,
142  const double *f, int L, double tau, int N)
143 {
144  assert(L > 0);
145  assert(N > 1);
146  int verbosity = 0;
147  int n, offset_lm, offset_r;
148  int flmsize = ssht_flm_size(L);
149  int frsize = ssht_fr_size_mw(L);
150  ssht_dl_method_t dl_method = SSHT_DL_TRAPANI;
151 
152  complex double *flmr;
153  flag_core_allocate_flmn(&flmr, L, N);
154 
155  for (n = 0; n < N; n++){
156  //printf("> Analysis: layer %i on %i\n",n+1,N);
157  offset_lm = n * flmsize;
158  offset_r = n * frsize;
159  ssht_core_mw_forward_sov_conv_sym_real(flmr + offset_lm, f + offset_r, L, dl_method, verbosity);
160  }
161 
162  double *nodes = (double*)calloc(N, sizeof(double));
163  double *weights = (double*)calloc(N, sizeof(double));
164  assert(nodes != NULL);
165  assert(weights != NULL);
166  flag_spherlaguerre_sampling(nodes, weights, tau, N);
167  //printf("> Mapped spherical Laguerre transform...");
168  fflush(NULL);
169  flag_spherlaguerre_mapped_analysis(flmn, flmr, nodes, weights, tau, N, flmsize);
170  //printf("done\n");
171  free(nodes);
172  free(weights);
173  free(flmr);
174 }
175 
176 void flag_core_synthesis_real(double *f,
177  const complex double *flmn,
178  const double *nodes, int Nnodes,
179  int L, double tau, int N)
180 {
181  assert(L > 0);
182  assert(N > 1);
183  int verbosity = 0;
184  int n, offset_lm, offset_r;
185  int flmsize = ssht_flm_size(L);
186  int frsize = ssht_fr_size_mw(L);
187  ssht_dl_method_t dl_method = SSHT_DL_TRAPANI;
188 
189  complex double *flmr;
190  //printf("> Mapped spherical Laguerre transform...");fflush(NULL);
191  flag_core_allocate_flmn(&flmr, L, Nnodes);
192  flag_spherlaguerre_mapped_synthesis(flmr, flmn, nodes, Nnodes, tau, N, flmsize);
193  //printf("done\n");
194 
195  for (n = 0; n < Nnodes; n++){
196  //printf("> Synthesis: layer %i on %i\n",n+1,N);
197  offset_lm = n * flmsize;
198  offset_r = n * frsize;
199  ssht_core_mw_inverse_sov_sym_real(f + offset_r, flmr + offset_lm, L, dl_method, verbosity);
200  }
201 
202  free(flmr);
203 }
204 
205 
206 double j_ell(double X, int l)
207 {
208  int L = l;
209  double JL, AX, AX2;
210  double LN2 = 0.6931471805599453094;
211  double ONEMLN2 = 0.30685281944005469058277;
212  double PID2 = 1.5707963267948966192313217;
213  double PID4 = 0.78539816339744830961566084582;
214  double ROOTPI12 = 21.269446210866192327578;
215  double GAMMA1 = 2.6789385347077476336556;
216  double GAMMA2 = 1.3541179394264004169452;
217  //double PI = 3.141592653589793238463;
218  double NU, NU2, BETA, BETA2, COSB;
219  double SX, SX2;
220  double COTB, COT3B, COT6B, SECB, SEC2B;
221  double TRIGARG, EXPTERM, L3;
222 
223  AX = pow(pow(X, 2.0), 0.5);
224  AX2 = pow(X, 2.0);
225 
226  switch(l){
227  case 0 :
228  if( AX < 0.1 )
229  JL = 1.0 - AX2 / 6.0 * (1.0 - AX2 / 20.0);
230  else
231  JL = sin(AX) / AX;
232  break;
233  case 1 :
234  if( AX < 0.2 )
235  JL = AX / 3.0 * (1.0 - AX2 / 10.0 * ( 1.0 - AX2 / 28.0 ) );
236  else
237  JL = (sin(AX) / AX - cos(AX)) / AX;
238  break;
239  case 2 :
240  if( AX < 0.3 )
241  JL = AX2 / 15.0 * (1.0 - AX2 / 14.0 * (1.0 - AX2 / 36.0));
242  else
243  JL=(-3.00 * cos(AX) / AX - sin(AX) * (1.0 - 3.0 / AX2)) / AX;
244  break;
245  case 3 :
246  if( AX < 0.4 )
247  JL = AX * AX2 / 105.0 * (1.0 - AX2 / 18.0 * (1.0 - AX2 / 44.0));
248  else
249  JL = (cos(AX) * (1.0 - 15.0 / AX2) - sin(AX) * (6.0 - 15.0 / AX2) / AX) / AX;
250  break;
251  case 4 :
252  if( AX < 0.6 )
253  JL = pow(AX2 ,2.0) / 945.0 * (1.0 - AX2 / 22.0 * (1.0 - AX2 / 52.0));
254  else
255  JL = (sin(AX) * (1.0 - (45.0 - 105.0 / AX2) / AX2) + cos(AX) * (10.0 - 105.0 / AX2) / AX) / AX;
256  break;
257  case 5 :
258  if( AX < 1.0 )
259  JL = pow(AX2, 2.0) * AX / 10395.0 * (1.0 - AX2 / 26.0 * (1.0 - AX2 / 60.0));
260  else
261  JL = (sin(AX) * (15.0 - (420.0 - 945.0 / AX2) / AX2) / AX - cos(AX) * (1.0 - (105.0 - 945.0 / AX2) / AX2)) / AX;
262  break;
263  case 6 :
264  if( AX < 1.0 )
265  JL = pow(AX2, 3.0) / 135135.0 * (1.0 - AX2 / 30.0 * (1.0 - AX2 / 68.0));
266  else
267  JL = (sin(AX) * ( - 1.0 + (210.0 - (4725.0 - 10395.0 / AX2) / AX2) / AX2) +
268  cos(AX) * ( - 21.0 + (1260.0 - 10395.0 / AX2) / AX2) / AX) / AX;
269  break;
270 
271  default :
272 
273  NU = 0.50 + L;
274  NU2 = pow(NU, 2.0);
275  if(AX < 1.0 - 40)
276  JL = 0.0;
277  else if( (AX2 / (double)L) < 0.5){
278  JL = exp(L * log(AX / NU) - LN2 + NU * ONEMLN2 - (1.0 - (1.0 - 3.50 / NU2) / NU2 / 30.0) / 12.0 / NU)
279  / NU * (1.0 - AX2 / (4.0 * NU + 4.0) * (1.0 - AX2 / (8.0 * NU + 16.0) * (1.0 - AX2 / (12.0 * NU + 36.0))));
280  }else if((pow(L, 2.0) / AX) < 5.0 - 1){
281  BETA = AX - PID2 * (L + 1);
282  JL = (cos(BETA) * (1.0 - (NU2 - 0.250) * (NU2 - 2.250) / 8.0 / AX2 * (1.0 - (NU2 - 6.25) * (NU2 - 12.250) / 48.0 / AX2))
283  - sin(BETA) * (NU2 - 0.250) / 2.0 / AX * (1.0 - (NU2 - 2.250) * (NU2 - 6.250) / 24.0 / AX2 * (1.0 - (NU2 - 12.25) *
284  (NU2 - 20.25) / 80.0 / AX2)) ) / AX ;
285  }else{
286  L3 = pow(NU, 0.325);
287  if(AX < NU - 1.31 * L3){
288  COSB = NU / AX;
289  SX = sqrt(NU2 - AX2);
290  COTB = NU / SX;
291  SECB = AX / NU;
292  BETA = log(COSB + SX / AX);
293  COT3B = pow(COTB, 3.0);
294  COT6B = pow(COT3B, 2.0);
295  SEC2B = pow(SECB, 2.0);
296  EXPTERM = ( (2.0 + 3.0 * SEC2B) * COT3B / 24.0
297  - ( (4.0 + SEC2B) * SEC2B * COT6B / 16.0
298  + ((16.0 - (1512.0 + (3654.0 + 375.0 * SEC2B) * SEC2B) * SEC2B) * COT3B / 5760.0
299  + (32.0 + (288.0 + (232.0 + 13.0 * SEC2B) * SEC2B) * SEC2B) * SEC2B * COT6B / 128.0 / NU) * COT6B / NU)
300  / NU) / NU;
301  JL = sqrt(COTB * COSB) / (2.0 * NU) * exp( - NU * BETA + NU / COTB - EXPTERM);
302  }else if (AX > NU + 1.48 * L3){
303  COSB = NU / AX;
304  SX = sqrt(AX2 - NU2);
305  COTB = NU / SX;
306  SECB = AX / NU;
307  BETA = acos(COSB);
308  COT3B = pow(COTB, 3.0);
309  COT6B = pow(COT3B, 2.0);
310  SEC2B = pow(SECB, 2.0);
311  TRIGARG = NU / COTB - NU * BETA - PID4
312  - ((2.0 + 3.0 * SEC2B) * COT3B / 24.0
313  + (16.0 - (1512.0 + (3654.0 + 375.0 * SEC2B) * SEC2B) * SEC2B) * COT3B * COT6B / 5760.0 / NU2) / NU;
314  EXPTERM = ( (4.0 + SEC2B) * SEC2B * COT6B / 16.0
315  - (32.0 + (288.0 + (232.0 + 13.0 * SEC2B) * SEC2B) * SEC2B) * SEC2B * pow(COT6B, 2.0) / 128.0 / NU2) / NU2;
316  JL = sqrt(COTB * COSB) / NU * exp( - EXPTERM) * cos(TRIGARG);
317  }else{
318  BETA = AX - NU;
319  BETA2 = pow(BETA, 2.0);
320  SX = 6.0 / AX;
321  SX2 = pow(SX, 2.0);
322  SECB = pow(SX, 0.3333333333333333);
323  SEC2B = pow(SECB, 2.0);
324  JL = ( GAMMA1 * SECB + BETA * GAMMA2 * SEC2B
325  - (BETA2 / 18.0 - 1.0 / 45.0) * BETA * SX * SECB * GAMMA1
326  - ((BETA2 - 1.0) * BETA2 / 36.0 + 1.0 / 420.0) * SX * SEC2B * GAMMA2
327  + (((BETA2 / 1620.0 - 7.0 / 3240.0) * BETA2 + 1.0 / 648.0) * BETA2 - 1.0 / 8100.0) * SX2 * SECB * GAMMA1
328  + (((BETA2 / 4536.0 - 1.0 / 810.0) * BETA2 + 19.0 / 11340.0) * BETA2 - 13.0 / 28350.0) * BETA * SX2 * SEC2B * GAMMA2
329  - ((((BETA2 / 349920.0 - 1.0 / 29160.0) * BETA2 + 71.0 / 583200.0) * BETA2 - 121.0 / 874800.0) *
330  BETA2 + 7939.0 / 224532000.0) * BETA * SX2 * SX * SECB * GAMMA1) * sqrt(SX) / ROOTPI12;
331  }
332  }
333  break;
334  }
335 
336  if( (X < 0) && (L % 2 != 0) )
337  JL = -JL;
338 
339  return(JL);
340 }
341 
342 
343 void flag_spherbessel_basis(double *jell, const int ell, const double *nodes, int Nnodes)
344 {
345  assert(Nnodes > 0);
346  int i;
347  for (i = 0; i < Nnodes; i++)
348  jell[i] = j_ell(nodes[i], ell);
349 }
350 
351 
352 
353 
354 void flag_core_fourierbessel_analysis(complex double *flmn,
355  const complex double *f,
356  int L, double tau, int N)
357 {
358  assert(L > 0);
359  assert(N > 1);
360  //const int alpha = ALPHA;
361  int spin = 0;
362  int verbosity = 0;
363  int n;
364  int flmsize = ssht_flm_size(L);
365  int frsize = ssht_fr_size_mw(L);
366  ssht_dl_method_t dl_method = SSHT_DL_TRAPANI;
367  int offset_lm, offset_r;
368 
369  complex double *flmr;
370  flag_core_allocate_flmn(&flmr, L, N);
371 
372  for (n = 0; n < N; n++){
373  //printf("> Analysis: layer %i on %i\n", n+1,N+1);
374  offset_lm = n * flmsize;
375  offset_r = n * frsize;
376  ssht_core_mw_forward_sov_conv_sym(flmr + offset_lm, f + offset_r, L, spin, dl_method, verbosity);
377  }
378 
379  double *nodes = (double*)calloc(N, sizeof(double));
380  double *weights = (double*)calloc(N, sizeof(double));
381  assert(nodes != NULL);
382  assert(weights != NULL);
383  flag_spherlaguerre_sampling(nodes, weights, tau, N);
384  //printf("> Mapped spherical Laguerre transform...");
385  fflush(NULL);
386  flag_spherlaguerre_mapped_analysis(flmn, flmr, nodes, weights, tau, N, flmsize);
387  //printf("done\n");
388  free(nodes);
389  free(weights);
390  free(flmr);
391 
392 }
393 
394 void flag_core_fourierbessel_synthesis(complex double *f,
395  const complex double *flmn,
396  const double *nodes, int Nnodes,
397  int L, double tau, int N)
398 {
399  assert(L > 0);
400  assert(N > 1);
401  assert(nodes != NULL);
402  //const int alpha = ALPHA;
403  int spin = 0;
404  int verbosity = 0;
405  int n, offset_lm, offset_r;
406  int flmsize = ssht_flm_size(L);
407  int frsize = ssht_fr_size_mw(L);
408  ssht_dl_method_t dl_method = SSHT_DL_TRAPANI;
409 
410  complex double *flmr;
411  flag_core_allocate_flmn(&flmr, L, Nnodes);
412  //printf("> Mapped spherical Laguerre transform...");fflush(NULL);
413  flag_spherlaguerre_mapped_synthesis(flmr, flmn, nodes, Nnodes, tau, N, flmsize);
414  //printf("done\n");
415 
416  for (n = 0; n < Nnodes; n++){
417  //printf("> Synthesis: layer %i on %i\n",n+1,Nnodes);
418  offset_lm = n * flmsize;
419  offset_r = n * frsize;
420  ssht_core_mw_inverse_sov_sym(f + offset_r, flmr + offset_lm, L, spin, dl_method, verbosity);
421  }
422 
423  free(flmr);
424 }
flag_core_synthesis
void flag_core_synthesis(complex double *f, const complex double *flmn, const double *nodes, int Nnodes, int L, double tau, int N, int spin)
Definition: flag_core.c:110
flag_core_allocate_f
void flag_core_allocate_f(complex double **f, int L, int N)
Definition: flag_core.c:51
flag_core_allocate_flmn
void flag_core_allocate_flmn(complex double **flmn, int L, int N)
Definition: flag_core.c:41
flag_core_analysis
void flag_core_analysis(complex double *flmn, const complex double *f, int L, double tau, int N, int spin)
Definition: flag_core.c:71
flag_core_flmn_size
int flag_core_flmn_size(int L, int N)
Definition: flag_core.c:27
flag_core_synthesis_real
void flag_core_synthesis_real(double *f, const complex double *flmn, const double *nodes, int Nnodes, int L, double tau, int N)
Definition: flag_core.c:176
flag_spherlaguerre_mapped_synthesis
void flag_spherlaguerre_mapped_synthesis(complex double *f, const complex double *fn, const double *nodes, int Nnodes, double tau, int N, int mapsize)
Definition: flag_spherlaguerre.c:389
flag_core_fourierbessel_analysis
void flag_core_fourierbessel_analysis(complex double *flmn, const complex double *f, int L, double tau, int N)
Definition: flag_core.c:354
ssht_flm_size
int ssht_flm_size(int L)
Definition: flag_core.c:20
flag_core_allocate_f_real
void flag_core_allocate_f_real(double **f, int L, int N)
Definition: flag_core.c:61
flag_spherbessel_basis
void flag_spherbessel_basis(double *jell, const int ell, const double *nodes, int Nnodes)
Definition: flag_core.c:343
flag_spherlaguerre_mapped_analysis
void flag_spherlaguerre_mapped_analysis(complex double *fn, const complex double *f, const double *nodes, const double *weights, double tau, int N, int mapsize)
Definition: flag_spherlaguerre.c:337
ssht_fr_size_mw
int ssht_fr_size_mw(int L)
Definition: flag_core.c:13
j_ell
double j_ell(double X, int l)
Definition: flag_core.c:206
flag_core_fourierbessel_synthesis
void flag_core_fourierbessel_synthesis(complex double *f, const complex double *flmn, const double *nodes, int Nnodes, int L, double tau, int N)
Definition: flag_core.c:394
flag_spherlaguerre_sampling
void flag_spherlaguerre_sampling(double *nodes, double *weights, double tau, int N)
Definition: flag_spherlaguerre.c:236
flag_core_analysis_real
void flag_core_analysis_real(complex double *flmn, const double *f, int L, double tau, int N)
Definition: flag_core.c:141
flag_core_f_size_mw
int flag_core_f_size_mw(int L, int N)
Definition: flag_core.c:34
Generated on Tue Sep 22 2015 17:01:04 for FLAG by   doxygen 1.8.10