Skip to content

Commit 123c863

Browse files
ixkaelixkael
committed
Various updates to Flag and bessel features and tests, including alpha parameter set to zero
1 parent f6b5c00 commit 123c863

15 files changed

+820
-224
lines changed

‎include/flag.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,8 @@
11
#ifndef FLAG
22
#define FLAG
33

4+
#define ALPHA 0
5+
46
#include "flag_core.h"
57
#include "flag_sampling.h"
68
#include "flag_spherlaguerre.h"

‎include/flag_core.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77

88
#include <complex.h>
99

10+
1011
/*!
1112
* Allocate FLAG coefficients.
1213
*
@@ -77,6 +78,10 @@ void flag_core_analysis(complex double *flmn, const complex double *f, double R,
7778
*/
7879
void flag_core_synthesis(complex double *f, const complex double *flmn, const double *nodes, int Nnodes, int L, int N);
7980

81+
void flag_core_fourierbessel_analysis(complex double *flmn, const complex double *f, double R, int L, int N);
82+
void flag_core_fourierbessel_synthesis(complex double *f, const complex double *flmn, const double *nodes, int Nnodes, int L, int N);
83+
84+
8085
/*!
8186
* Perform Fourier-Laguerre analysis (MW sampling, real signal).
8287
*

‎include/flag_spherlaguerre.h‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22
#ifndef FLAG_SPHERLAGUERRE
33
#define FLAG_SPHERLAGUERRE
44

5-
#include <complex.h>
5+
#include <complex.h>
66

77
/*!
88
* Compute Gauss-Laguerre quadrature (nodes and weights).
@@ -34,6 +34,8 @@ double flag_spherlaguerre_tau(double R, int N);
3434
*/
3535
void flag_spherlaguerre_sampling(double *nodes, double *weights, double R, int N);
3636

37+
void flag_spherbessel_sampling(double *nodes, double *weights, double R, int N);
38+
3739
/*!
3840
* Perform spherical Laguerre analysis.
3941
*
@@ -57,6 +59,9 @@ void flag_spherlaguerre_analysis(double *fn, const double *f, const double *node
5759
*/
5860
void flag_spherlaguerre_synthesis(double *f, const double *fn, const double *nodes, int Nnodes, int N);
5961

62+
void flag_spherlaguerre_synthesis_gen(double *f, const double *fn, const double *nodes, int Nnodes, int N, int alpha);
63+
64+
6065
/*!
6166
* Allocate spherical Laguerre sampling scheme.
6267
*

‎src/main/c/flag_core.c‎

Lines changed: 114 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,12 @@
11
// FLAG package
2-
// Copyright (C) 2012
2+
// Copyright (C) 2012
33
// Boris Leistedt & Jason McEwen
44

55
#include "flag.h"
66
#include <math.h>
77
#include <stdlib.h>
8-
#include <complex.h>
9-
#include <fftw3.h>
8+
#include <complex.h>
9+
#include <fftw3.h>
1010
#include <ssht.h>
1111
#include <assert.h>
1212

@@ -68,31 +68,31 @@ void flag_core_allocate_f_real(double **f, int L, int N)
6868
assert(f != NULL);
6969
}
7070

71-
void flag_core_analysis(complex double *flmn,
72-
const complex double *f,
71+
void flag_core_analysis(complex double *flmn,
72+
const complex double *f,
7373
double R, int L, int N)
7474
{
7575
assert(L > 0);
7676
assert(N > 1);
77-
//const int alpha = 2;
77+
//const int alpha = ALPHA;
7878
int spin = 0;
7979
int verbosity = 0;
8080
int n;
8181
int flmsize = ssht_flm_size(L);
8282
int frsize = ssht_fr_size_mw(L);
8383
ssht_dl_method_t dl_method = SSHT_DL_TRAPANI;
8484
int offset_lm, offset_r;
85-
85+
8686
complex double *flmr;
8787
flag_core_allocate_flmn(&flmr, L, N);
88-
88+
8989
for (n = 0; n < N; n++){
9090
//printf("> Analysis: layer %i on %i\n", n+1,N+1);
9191
offset_lm = n * flmsize;
9292
offset_r = n * frsize;
9393
ssht_core_mw_forward_sov_conv_sym(flmr + offset_lm, f + offset_r, L, spin, dl_method, verbosity);
9494
}
95-
95+
9696
double *nodes = (double*)calloc(N, sizeof(double));
9797
double *weights = (double*)calloc(N, sizeof(double));
9898
assert(nodes != NULL);
@@ -108,15 +108,15 @@ void flag_core_analysis(complex double *flmn,
108108

109109
}
110110

111-
void flag_core_synthesis(complex double *f,
112-
const complex double *flmn,
111+
void flag_core_synthesis(complex double *f,
112+
const complex double *flmn,
113113
const double *nodes, int Nnodes,
114114
int L, int N)
115115
{
116116
assert(L > 0);
117117
assert(N > 1);
118118
assert(nodes != NULL);
119-
//const int alpha = 2;
119+
//const int alpha = ALPHA;
120120
int spin = 0;
121121
int verbosity = 0;
122122
int n, offset_lm, offset_r;
@@ -129,7 +129,7 @@ void flag_core_synthesis(complex double *f,
129129
//printf("> Mapped spherical Laguerre transform...");fflush(NULL);
130130
flag_spherlaguerre_mapped_synthesis(flmr, flmn, nodes, Nnodes, N, flmsize);
131131
//printf("done\n");
132-
132+
133133
for (n = 0; n < Nnodes; n++){
134134
//printf("> Synthesis: layer %i on %i\n",n+1,Nnodes);
135135
offset_lm = n * flmsize;
@@ -140,7 +140,7 @@ void flag_core_synthesis(complex double *f,
140140
free(flmr);
141141
}
142142

143-
void flag_core_analysis_real(complex double *flmn,
143+
void flag_core_analysis_real(complex double *flmn,
144144
const double *f, double R, int L, int N)
145145
{
146146
assert(L > 0);
@@ -175,9 +175,9 @@ void flag_core_analysis_real(complex double *flmn,
175175
free(flmr);
176176
}
177177

178-
void flag_core_synthesis_real(double *f,
179-
const complex double *flmn,
180-
const double *nodes, int Nnodes,
178+
void flag_core_synthesis_real(double *f,
179+
const complex double *flmn,
180+
const double *nodes, int Nnodes,
181181
int L, int N)
182182
{
183183
assert(L > 0);
@@ -232,13 +232,13 @@ double j_ell(double X, int l)
232232
else
233233
JL = sin(AX) / AX;
234234
break;
235-
case 1 :
235+
case 1 :
236236
if( AX < 0.2 )
237237
JL = AX / 3.0 * (1.0 - AX2 / 10.0 * ( 1.0 - AX2 / 28.0 ) );
238238
else
239239
JL = (sin(AX) / AX - cos(AX)) / AX;
240240
break;
241-
case 2 :
241+
case 2 :
242242
if( AX < 0.3 )
243243
JL = AX2 / 15.0 * (1.0 - AX2 / 14.0 * (1.0 - AX2 / 36.0));
244244
else
@@ -266,23 +266,23 @@ double j_ell(double X, int l)
266266
if( AX < 1.0 )
267267
JL = pow(AX2, 3.0) / 135135.0 * (1.0 - AX2 / 30.0 * (1.0 - AX2 / 68.0));
268268
else
269-
JL = (sin(AX) * ( - 1.0 + (210.0 - (4725.0 - 10395.0 / AX2) / AX2) / AX2) +
269+
JL = (sin(AX) * ( - 1.0 + (210.0 - (4725.0 - 10395.0 / AX2) / AX2) / AX2) +
270270
cos(AX) * ( - 21.0 + (1260.0 - 10395.0 / AX2) / AX2) / AX) / AX;
271271
break;
272272

273273
default :
274-
274+
275275
NU = 0.50 + L;
276276
NU2 = pow(NU, 2.0);
277277
if(AX < 1.0 - 40)
278278
JL = 0.0;
279279
else if( (AX2 / (double)L) < 0.5){
280-
JL = exp(L * log(AX / NU) - LN2 + NU * ONEMLN2 - (1.0 - (1.0 - 3.50 / NU2) / NU2 / 30.0) / 12.0 / NU)
280+
JL = exp(L * log(AX / NU) - LN2 + NU * ONEMLN2 - (1.0 - (1.0 - 3.50 / NU2) / NU2 / 30.0) / 12.0 / NU)
281281
/ 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))));
282282
}else if((pow(L, 2.0) / AX) < 5.0 - 1){
283283
BETA = AX - PID2 * (L + 1);
284-
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))
285-
- 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+
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))
285+
- sin(BETA) * (NU2 - 0.250) / 2.0 / AX * (1.0 - (NU2 - 2.250) * (NU2 - 6.250) / 24.0 / AX2 * (1.0 - (NU2 - 12.25) *
286286
(NU2 - 20.25) / 80.0 / AX2)) ) / AX ;
287287
}else{
288288
L3 = pow(NU, 0.325);
@@ -295,10 +295,10 @@ double j_ell(double X, int l)
295295
COT3B = pow(COTB, 3.0);
296296
COT6B = pow(COT3B, 2.0);
297297
SEC2B = pow(SECB, 2.0);
298-
EXPTERM = ( (2.0 + 3.0 * SEC2B) * COT3B / 24.0
299-
- ( (4.0 + SEC2B) * SEC2B * COT6B / 16.0
300-
+ ((16.0 - (1512.0 + (3654.0 + 375.0 * SEC2B) * SEC2B) * SEC2B) * COT3B / 5760.0
301-
+ (32.0 + (288.0 + (232.0 + 13.0 * SEC2B) * SEC2B) * SEC2B) * SEC2B * COT6B / 128.0 / NU) * COT6B / NU)
298+
EXPTERM = ( (2.0 + 3.0 * SEC2B) * COT3B / 24.0
299+
- ( (4.0 + SEC2B) * SEC2B * COT6B / 16.0
300+
+ ((16.0 - (1512.0 + (3654.0 + 375.0 * SEC2B) * SEC2B) * SEC2B) * COT3B / 5760.0
301+
+ (32.0 + (288.0 + (232.0 + 13.0 * SEC2B) * SEC2B) * SEC2B) * SEC2B * COT6B / 128.0 / NU) * COT6B / NU)
302302
/ NU) / NU;
303303
JL = sqrt(COTB * COSB) / (2.0 * NU) * exp( - NU * BETA + NU / COTB - EXPTERM);
304304
}else if (AX > NU + 1.48 * L3){
@@ -310,10 +310,10 @@ double j_ell(double X, int l)
310310
COT3B = pow(COTB, 3.0);
311311
COT6B = pow(COT3B, 2.0);
312312
SEC2B = pow(SECB, 2.0);
313-
TRIGARG = NU / COTB - NU * BETA - PID4
314-
- ((2.0 + 3.0 * SEC2B) * COT3B / 24.0
313+
TRIGARG = NU / COTB - NU * BETA - PID4
314+
- ((2.0 + 3.0 * SEC2B) * COT3B / 24.0
315315
+ (16.0 - (1512.0 + (3654.0 + 375.0 * SEC2B) * SEC2B) * SEC2B) * COT3B * COT6B / 5760.0 / NU2) / NU;
316-
EXPTERM = ( (4.0 + SEC2B) * SEC2B * COT6B / 16.0
316+
EXPTERM = ( (4.0 + SEC2B) * SEC2B * COT6B / 16.0
317317
- (32.0 + (288.0 + (232.0 + 13.0 * SEC2B) * SEC2B) * SEC2B) * SEC2B * pow(COT6B, 2.0) / 128.0 / NU2) / NU2;
318318
JL = sqrt(COTB * COSB) / NU * exp( - EXPTERM) * cos(TRIGARG);
319319
}else{
@@ -323,19 +323,19 @@ double j_ell(double X, int l)
323323
SX2 = pow(SX, 2.0);
324324
SECB = pow(SX, 0.3333333333333333);
325325
SEC2B = pow(SECB, 2.0);
326-
JL = ( GAMMA1 * SECB + BETA * GAMMA2 * SEC2B
327-
- (BETA2 / 18.0 - 1.0 / 45.0) * BETA * SX * SECB * GAMMA1
328-
- ((BETA2 - 1.0) * BETA2 / 36.0 + 1.0 / 420.0) * SX * SEC2B * GAMMA2
329-
+ (((BETA2 / 1620.0 - 7.0 / 3240.0) * BETA2 + 1.0 / 648.0) * BETA2 - 1.0 / 8100.0) * SX2 * SECB * GAMMA1
330-
+ (((BETA2 / 4536.0 - 1.0 / 810.0) * BETA2 + 19.0 / 11340.0) * BETA2 - 13.0 / 28350.0) * BETA * SX2 * SEC2B * GAMMA2
331-
- ((((BETA2 / 349920.0 - 1.0 / 29160.0) * BETA2 + 71.0 / 583200.0) * BETA2 - 121.0 / 874800.0) *
326+
JL = ( GAMMA1 * SECB + BETA * GAMMA2 * SEC2B
327+
- (BETA2 / 18.0 - 1.0 / 45.0) * BETA * SX * SECB * GAMMA1
328+
- ((BETA2 - 1.0) * BETA2 / 36.0 + 1.0 / 420.0) * SX * SEC2B * GAMMA2
329+
+ (((BETA2 / 1620.0 - 7.0 / 3240.0) * BETA2 + 1.0 / 648.0) * BETA2 - 1.0 / 8100.0) * SX2 * SECB * GAMMA1
330+
+ (((BETA2 / 4536.0 - 1.0 / 810.0) * BETA2 + 19.0 / 11340.0) * BETA2 - 13.0 / 28350.0) * BETA * SX2 * SEC2B * GAMMA2
331+
- ((((BETA2 / 349920.0 - 1.0 / 29160.0) * BETA2 + 71.0 / 583200.0) * BETA2 - 121.0 / 874800.0) *
332332
BETA2 + 7939.0 / 224532000.0) * BETA * SX2 * SX * SECB * GAMMA1) * sqrt(SX) / ROOTPI12;
333333
}
334334
}
335335
break;
336336
}
337-
338-
if( (X < 0) && (L % 2 != 0) )
337+
338+
if( (X < 0) && (L % 2 != 0) )
339339
JL = -JL;
340340

341341
return(JL);
@@ -349,3 +349,78 @@ void flag_spherbessel_basis(double *jell, const int ell, const double *nodes, in
349349
for (i = 0; i < Nnodes; i++)
350350
jell[i] = j_ell(nodes[i], ell);
351351
}
352+
353+
354+
355+
356+
void flag_core_fourierbessel_analysis(complex double *flmn,
357+
const complex double *f,
358+
double R, int L, int N)
359+
{
360+
assert(L > 0);
361+
assert(N > 1);
362+
//const int alpha = ALPHA;
363+
int spin = 0;
364+
int verbosity = 0;
365+
int n;
366+
int flmsize = ssht_flm_size(L);
367+
int frsize = ssht_fr_size_mw(L);
368+
ssht_dl_method_t dl_method = SSHT_DL_TRAPANI;
369+
int offset_lm, offset_r;
370+
371+
complex double *flmr;
372+
flag_core_allocate_flmn(&flmr, L, N);
373+
374+
for (n = 0; n < N; n++){
375+
//printf("> Analysis: layer %i on %i\n", n+1,N+1);
376+
offset_lm = n * flmsize;
377+
offset_r = n * frsize;
378+
ssht_core_mw_forward_sov_conv_sym(flmr + offset_lm, f + offset_r, L, spin, dl_method, verbosity);
379+
}
380+
381+
double *nodes = (double*)calloc(N, sizeof(double));
382+
double *weights = (double*)calloc(N, sizeof(double));
383+
assert(nodes != NULL);
384+
assert(weights != NULL);
385+
flag_spherlaguerre_sampling(nodes, weights, R, N);
386+
//printf("> Mapped spherical Laguerre transform...");
387+
fflush(NULL);
388+
flag_spherlaguerre_mapped_analysis(flmn, flmr, nodes, weights, N, flmsize);
389+
//printf("done\n");
390+
free(nodes);
391+
free(weights);
392+
free(flmr);
393+
394+
}
395+
396+
void flag_core_fourierbessel_synthesis(complex double *f,
397+
const complex double *flmn,
398+
const double *nodes, int Nnodes,
399+
int L, int N)
400+
{
401+
assert(L > 0);
402+
assert(N > 1);
403+
assert(nodes != NULL);
404+
//const int alpha = ALPHA;
405+
int spin = 0;
406+
int verbosity = 0;
407+
int n, offset_lm, offset_r;
408+
int flmsize = ssht_flm_size(L);
409+
int frsize = ssht_fr_size_mw(L);
410+
ssht_dl_method_t dl_method = SSHT_DL_TRAPANI;
411+
412+
complex double *flmr;
413+
flag_core_allocate_flmn(&flmr, L, Nnodes);
414+
//printf("> Mapped spherical Laguerre transform...");fflush(NULL);
415+
flag_spherlaguerre_mapped_synthesis(flmr, flmn, nodes, Nnodes, N, flmsize);
416+
//printf("done\n");
417+
418+
for (n = 0; n < Nnodes; n++){
419+
//printf("> Synthesis: layer %i on %i\n",n+1,Nnodes);
420+
offset_lm = n * flmsize;
421+
offset_r = n * frsize;
422+
ssht_core_mw_inverse_sov_sym(f + offset_r, flmr + offset_lm, L, spin, dl_method, verbosity);
423+
}
424+
425+
free(flmr);
426+
}

0 commit comments

Comments
 (0)