FLAGLET: src/main/c/flaglet_axisym.c Source File

Go to the documentation of this file.
00001 // FLAGLET package
00002 // Copyright (C) 2012 
00003 // Boris Leistedt & Jason McEwen
00004 
00005 #include "flaglet.h"
00006 #include <stdlib.h>
00007 #include <math.h>
00008 #include <complex.h> 
00009 #include <flag.h>
00010 #include <s2let.h>
00011 #include <assert.h>
00012 
00024 void flaglet_axisym_allocate_wav_lmp(double **wav_lmp, double **scal_lmp, int B_l, int B_p, int L, int P)
00025 {
00026     int J_l = s2let_j_max(L, B_l);
00027     int J_p = s2let_j_max(P, B_p);
00028     *wav_lmp = (double*)calloc( (J_l+1) * L * L * (J_p+1) * P, sizeof(double));
00029     *scal_lmp = (double*)calloc( L * L * P, sizeof(double));
00030 }
00031 
00045 void flaglet_axisym_wav_lmp(double *wav_lmp, double *scal_lmp, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00046 {
00047     int jl, jp, l, n, m, indjjlmp, indlmp;
00048     int J_l = s2let_j_max(L, B_l);
00049     int J_p = s2let_j_max(P, B_p);
00050 
00051     double wav0, scal0;
00052 
00053     double *kappa_ln, *kappa0_ln;
00054     flaglet_tiling_axisym_allocate(&kappa_ln, &kappa0_ln, B_l, B_p, L, P);
00055     
00056     flaglet_tiling_axisym(kappa_ln, kappa0_ln, B_l, B_p, L, P, J_min_l, J_min_p);
00057 
00058     for (jp = J_min_p; jp <= J_p; jp++){
00059         for (jl = J_min_l; jl <= J_l; jl++){
00060             for (n = 0; n < P; n++){
00061                 for (l = 0; l < L; l++){
00062                     wav0 = sqrt( (2 * l + 1) / (4.0 * PI) ) *
00063                         kappa_ln[ jp*(J_l+1)*L*P  + jl*L*P + n*L + l ];
00064                     for (m = -l; m <= l ; m++){
00065                         indjjlmp = jjlmp2ind(jl,jp,l,m,n,J_l,J_p,L,P);
00066                         wav_lmp[indjjlmp] =  wav0;
00067                     }
00068                 }
00069             }
00070         }
00071     }
00072     for (n = 0; n < P; n++){
00073         for (l = 0; l < L; l++){
00074             scal0 = sqrt( (2 * l + 1) / (4.0 * PI) ) *
00075                 kappa0_ln[ n*L + l ];
00076             for (m = -l; m <= l ; m++){
00077                 indlmp = lmp2ind(l,m,n,L);
00078                 scal_lmp[indlmp] = scal0;
00079             }
00080         }
00081     }
00082 
00083     free(kappa_ln);
00084     free(kappa0_ln);
00085 }
00086 
00100 void flaglet_axisym_allocate_f_wav_lmp(complex double **f_wav_lmp, complex double **f_scal_lmp, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00101 {
00102     int J_l = s2let_j_max(L, B_l);
00103     int J_p = s2let_j_max(P, B_p);
00104     *f_wav_lmp = (complex double*)calloc( (J_l+1-J_min_l) * L * L * (J_p+1-J_min_p) * P, sizeof(complex double));
00105     *f_scal_lmp = (complex double*)calloc( L * L * P, sizeof(complex double));
00106 }
00107 
00121 void flaglet_axisym_allocate_f_wav(complex double **f_wav, complex double **f_scal, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00122 {
00123     int J_l = s2let_j_max(L, B_l);
00124     int J_p = s2let_j_max(P, B_p);
00125     *f_wav = (complex double*)calloc( (J_l+1-J_min_l) * L * (2*L-1) * (J_p+1-J_min_p) * P, sizeof(complex double));
00126     *f_scal = (complex double*)calloc( L * (2*L-1) * P, sizeof(complex double));
00127 }
00128 
00142 void flaglet_axisym_allocate_f_wav_real(double **f_wav, double **f_scal, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00143 {
00144     int J_l = s2let_j_max(L, B_l);
00145     int J_p = s2let_j_max(P, B_p);
00146     *f_wav = (double*)calloc( (J_l+1-J_min_l) * L * (2*L-1) * (J_p+1-J_min_p) * P, sizeof(double));
00147     *f_scal = (double*)calloc( L * (2*L-1) * P, sizeof(double));
00148 }
00149 
00163 void flaglet_axisym_allocate_f_wav_multires_lmp(complex double **f_wav_lmp, complex double **f_scal_lmp, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00164 {
00165     int J_l = s2let_j_max(L, B_l);
00166     int J_p = s2let_j_max(P, B_p);
00167     int jp, jl, bandlimit_p, bandlimit_l, total = 0;
00168     for (jp = J_min_p; jp <= J_p; jp++){
00169         bandlimit_p = MIN(s2let_bandlimit(B_p, jp), P);
00170         for (jl = J_min_l; jl <= J_l; jl++){
00171             bandlimit_l = MIN(s2let_bandlimit(B_l, jl), L);
00172             total += bandlimit_l * bandlimit_l * bandlimit_p;
00173         }
00174     }
00175     *f_wav_lmp = (complex double*)calloc( total, sizeof(complex double));
00176     *f_scal_lmp = (complex double*)calloc( L * L * P, sizeof(complex double));
00177 }
00178 
00192 void flaglet_axisym_allocate_f_wav_multires(complex double **f_wav, complex double **f_scal, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00193 {
00194     int J_l = s2let_j_max(L, B_l);
00195     int J_p = s2let_j_max(P, B_p);
00196     int jp, jl, bandlimit_p, bandlimit_l, total = 0;
00197     for (jp = J_min_p; jp <= J_p; jp++){
00198         bandlimit_p = MIN(s2let_bandlimit(B_p, jp), P);
00199         for (jl = J_min_l; jl <= J_l; jl++){
00200             bandlimit_l = MIN(s2let_bandlimit(B_l, jl), L);
00201             total += bandlimit_l * (2 * bandlimit_l - 1) * bandlimit_p;
00202         }
00203     }
00204     *f_wav = (complex double*)calloc( total, sizeof(complex double));
00205     *f_scal = (complex double*)calloc( L * (2*L-1) * P, sizeof(complex double));
00206 }
00207 
00221 void flaglet_axisym_allocate_f_wav_multires_real(double **f_wav, double **f_scal, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00222 {
00223     int J_l = s2let_j_max(L, B_l);
00224     int J_p = s2let_j_max(P, B_p);
00225     int jp, jl, bandlimit_p, bandlimit_l, total = 0;
00226     for (jp = J_min_p; jp <= J_p; jp++){
00227         bandlimit_p = MIN(s2let_bandlimit(B_p, jp), P);
00228         for (jl = J_min_l; jl <= J_l; jl++){
00229             bandlimit_l = MIN(s2let_bandlimit(B_l, jl), L);
00230             total += bandlimit_l * (2 * bandlimit_l - 1) * bandlimit_p;
00231         }
00232     }
00233     *f_wav = (double*)calloc( total, sizeof(double));
00234     *f_scal = (double*)calloc( L * (2*L-1) * P, sizeof(double));
00235 }
00236 
00254 void flaglet_axisym_wav_analysis_multires_lmp(complex double *f_wav_lmp, complex double *f_scal_lmp, const complex double *flmp, const double *wav_lmp, const double *scal_lmp, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00255 {
00256     int bandlimit_p, bandlimit_l, offset, jl, jp, l, m, n, indjjlmp, indlmp;
00257     int J_l = s2let_j_max(L, B_l);
00258     int J_p = s2let_j_max(P, B_p);
00259 
00260     offset = 0;
00261     for (jp = J_min_p; jp <= J_p; jp++){
00262         bandlimit_p = MIN(s2let_bandlimit(B_p, jp), P);
00263         for (jl = J_min_l; jl <= J_l; jl++){
00264             bandlimit_l = MIN(s2let_bandlimit(B_l, jl), L);
00265             for (n = 0; n < bandlimit_p; n++){
00266                 for (l = 0; l < bandlimit_l; l++){
00267                     for (m = -l; m <= l ; m++){
00268                         indjjlmp = jjlmp2ind(jl,jp,l,m,n,J_l,J_p,L,P);
00269                         indlmp = lmp2ind(l,m,n,L);
00270                         f_wav_lmp[offset + lmp2ind(l,m,n,bandlimit_l)] = flmp[indlmp] * sqrt((4.0*PI)/(2.0*l+1.0)) * wav_lmp[indjjlmp] ;
00271                     }
00272                 }
00273             }
00274             offset += bandlimit_l * bandlimit_l * bandlimit_p;
00275         }
00276     }
00277 
00278     for (n = 0; n < P; n++){
00279         for (l = 0; l < L; l++){
00280             for (m = -l; m <= l ; m++){
00281                 indlmp = lmp2ind(l,m,n,L);
00282                 f_scal_lmp[indlmp] = flmp[indlmp] * sqrt((4.0*PI)/(2.0*l+1.0)) * scal_lmp[indlmp] ;
00283             }
00284         }
00285     }
00286 }
00287 
00305 void flaglet_axisym_wav_synthesis_multires_lmp(complex double *flmp, const complex double *f_wav_lmp, const complex double *f_scal_lmp, const double *wav_lmp, const double *scal_lmp, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00306 {
00307     int bandlimit_p, bandlimit_l, offset, jl, jp, l, m, n, indjjlmp, indlmp;
00308     int J_l = s2let_j_max(L, B_l);
00309     int J_p = s2let_j_max(P, B_p);
00310 
00311     offset = 0;
00312     for (jp = J_min_p; jp <= J_p; jp++){
00313         bandlimit_p = MIN(s2let_bandlimit(B_p, jp), P);
00314         for (jl = J_min_l; jl <= J_l; jl++){
00315             bandlimit_l = MIN(s2let_bandlimit(B_l, jl), L);
00316             for (n = 0; n < bandlimit_p; n++){
00317                 for (l = 0; l < bandlimit_l; l++){
00318                     for (m = -l; m <= l ; m++){
00319                         indjjlmp = jjlmp2ind(jl,jp,l,m,n,J_l,J_p,L,P);
00320                         indlmp = lmp2ind(l,m,n,L);
00321                         flmp[indlmp] += f_wav_lmp[offset + lmp2ind(l,m,n,bandlimit_l)] * sqrt((4.0*PI)/(2.0*l+1.0)) * wav_lmp[indjjlmp] ;
00322                     }
00323                 }
00324             }
00325             offset += bandlimit_l * bandlimit_l * bandlimit_p;
00326         }
00327     }
00328 
00329     bandlimit_p = MIN(s2let_bandlimit(B_p, J_min_p-1), P);
00330     bandlimit_l = MIN(s2let_bandlimit(B_l, J_min_l-1), L);
00331     for (n = 0; n < P; n++){
00332         for (l = 0; l < L; l++){
00333             for (m = -l; m <= l ; m++){
00334                 indlmp = lmp2ind(l,m,n,L);
00335                 flmp[indlmp] += f_scal_lmp[indlmp] * sqrt((4.0*PI)/(2.0*l+1.0)) * scal_lmp[indlmp] ;
00336             }
00337         }
00338     }
00339     
00340 }
00341 
00359 void flaglet_axisym_wav_analysis_lmp(complex double *f_wav_lmp, complex double *f_scal_lmp, const complex double *flmp, const double *wav_lmp, const double *scal_lmp, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00360 {
00361     int offset, jl, jp, l, m, n, indjjlmp, indlmp;
00362     int J_l = s2let_j_max(L, B_l);
00363     int J_p = s2let_j_max(P, B_p);
00364 
00365     offset = 0;
00366     for (jp = J_min_p; jp <= J_p; jp++){
00367         for (jl = J_min_l; jl <= J_l; jl++){
00368             for (n = 0; n < P; n++){
00369                 for (l = 0; l < L; l++){
00370                     for (m = -l; m <= l ; m++){
00371                         indjjlmp = jjlmp2ind(jl,jp,l,m,n,J_l,J_p,L,P);
00372                         indlmp = lmp2ind(l,m,n,L);
00373                         f_wav_lmp[offset + indlmp] = flmp[indlmp] * sqrt((4.0*PI)/(2.0*l+1.0)) * wav_lmp[indjjlmp] ;
00374                     }
00375                 }
00376             }
00377             offset += L * L * P;
00378         }
00379     }
00380     for (n = 0; n < P; n++){
00381         for (l = 0; l < L; l++){
00382             for (m = -l; m <= l ; m++){
00383                 indlmp = lmp2ind(l,m,n,L);
00384                 f_scal_lmp[indlmp] = flmp[indlmp] * sqrt((4.0*PI)/(2.0*l+1.0)) * scal_lmp[indlmp] ;
00385             }
00386         }
00387     }
00388 }
00389 
00407 void flaglet_axisym_wav_synthesis_lmp(complex double *flmp, const complex double *f_wav_lmp, const complex double *f_scal_lmp, const double *wav_lmp, const double *scal_lmp, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00408 {
00409     int offset, jl, jp, l, m, n, indjjlmp, indlmp;
00410     int J_l = s2let_j_max(L, B_l);
00411     int J_p = s2let_j_max(P, B_p);
00412 
00413     offset = 0;
00414     for (jp = J_min_p; jp <= J_p; jp++){
00415         for (jl = J_min_l; jl <= J_l; jl++){
00416             for (n = 0; n < P; n++){
00417                 for (l = 0; l < L; l++){
00418                     for (m = -l; m <= l ; m++){
00419                         indjjlmp = jjlmp2ind(jl,jp,l,m,n,J_l,J_p,L,P);
00420                         indlmp = lmp2ind(l,m,n,L);
00421                         flmp[indlmp] += f_wav_lmp[offset + indlmp] * sqrt((4.0*PI)/(2.0*l+1.0)) * wav_lmp[indjjlmp] ;
00422                     }
00423                 }
00424             }
00425             offset += L * L * P;
00426         }
00427     }
00428 
00429     for (n = 0; n < P; n++){
00430         for (l = 0; l < L; l++){
00431             for (m = -l; m <= l ; m++){
00432                 indlmp = lmp2ind(l,m,n,L);
00433                 flmp[indlmp] += f_scal_lmp[indlmp] * sqrt((4.0*PI)/(2.0*l+1.0)) * scal_lmp[indlmp] ;
00434             }
00435         }
00436     }
00437     
00438 }
00439 
00440 
00459 void flaglet_axisym_wav_analysis_multires(complex double *f_wav, complex double *f_scal, const complex double *f, double R, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00460 {
00461     complex double *flmp;
00462     flag_core_allocate_flmn(&flmp, L, P);
00463     flag_core_analysis(flmp, f, R, L, P);
00464 
00465     double *wav_lmp, *scal_lmp;
00466     flaglet_axisym_allocate_wav_lmp(&wav_lmp, &scal_lmp, B_l, B_p, L, P);
00467     flaglet_axisym_wav_lmp(wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00468 
00469     complex double *f_wav_lmp, *f_scal_lmp;
00470     flaglet_axisym_allocate_f_wav_multires_lmp(&f_wav_lmp, &f_scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00471     flaglet_axisym_wav_analysis_multires_lmp(f_wav_lmp, f_scal_lmp, flmp, wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00472 
00473     double *nodes = (double*)calloc(P, sizeof(double));
00474     double *weights = (double*)calloc(P, sizeof(double));
00475     assert(nodes != NULL);
00476     assert(weights != NULL);
00477     flag_spherlaguerre_sampling(nodes, weights, R, P);
00478 
00479     int bandlimit_p, bandlimit_l, offset_lmp, offset, jl, jp;
00480     int J_l = s2let_j_max(L, B_l);
00481     int J_p = s2let_j_max(P, B_p);
00482 
00483     flag_core_synthesis(f_scal, f_scal_lmp, nodes, P, L, P);
00484 
00485     offset_lmp = 0;
00486     offset = 0;
00487     for (jp = J_min_p; jp <= J_p; jp++){
00488         bandlimit_p = MIN(s2let_bandlimit(B_p, jp), P);
00489         flag_spherlaguerre_sampling(nodes, weights, R, bandlimit_p);
00490         for (jl = J_min_l; jl <= J_l; jl++){
00491             bandlimit_l = MIN(s2let_bandlimit(B_l, jl), L);         
00492             //offset_lmp = jp * (J_l + 1) * L * L * P   +  jl * L * L * P;
00493             //offset = jp * (J_l + 1) * L * (2*L-1) * (P)  +  jl * L * (2*L-1) * (P);
00494             flag_core_synthesis(f_wav + offset, f_wav_lmp + offset_lmp, nodes, bandlimit_p, bandlimit_l, bandlimit_p);
00495             offset_lmp += bandlimit_l * bandlimit_l * bandlimit_p ;
00496             offset += bandlimit_l * (2 * bandlimit_l - 1) * bandlimit_p;
00497         }
00498     }
00499 
00500     free(nodes);
00501     free(weights);
00502     free(f_wav_lmp);
00503     free(f_scal_lmp);
00504     free(flmp);
00505 }
00506 
00525 void flaglet_axisym_wav_synthesis_multires(complex double *f, const complex double *f_wav, const complex double *f_scal, double R, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00526 {
00527     complex double *f_wav_lmp, *f_scal_lmp;
00528     flaglet_axisym_allocate_f_wav_multires_lmp(&f_wav_lmp, &f_scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00529 
00530     int bandlimit_p, bandlimit_l, offset_lmp, offset, jl, jp;
00531     int J_l = s2let_j_max(L, B_l);
00532     int J_p = s2let_j_max(P, B_p);
00533 
00534     flag_core_analysis(f_scal_lmp, f_scal, R, L, P);
00535 
00536     offset_lmp = 0;
00537     offset = 0;
00538     for (jp = J_min_p; jp <= J_p; jp++){
00539         bandlimit_p = MIN(s2let_bandlimit(B_p, jp), P);
00540         for (jl = J_min_l; jl <= J_l; jl++){
00541             bandlimit_l = MIN(s2let_bandlimit(B_l, jl), L);
00542             //offset_lmp = jp * (J_l + 1) * L * L * P   +  jl * L * L * P;
00543             //offset = jp * (J_l + 1) * L * (2*L-1) * (P)   +  jl * L * (2*L-1) * (P);
00544             flag_core_analysis(f_wav_lmp + offset_lmp, f_wav + offset, R, bandlimit_l, bandlimit_p);
00545             offset_lmp += bandlimit_l * bandlimit_l * bandlimit_p ;
00546             offset += bandlimit_l * (2 * bandlimit_l - 1) * bandlimit_p;
00547         }
00548     }
00549 
00550     double *wav_lmp, *scal_lmp ;
00551     flaglet_axisym_allocate_wav_lmp(&wav_lmp, &scal_lmp, B_l, B_p, L, P);
00552     flaglet_axisym_wav_lmp(wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00553 
00554     complex double *flmp;
00555     flag_core_allocate_flmn(&flmp, L, P);
00556     flaglet_axisym_wav_synthesis_multires_lmp(flmp, f_wav_lmp, f_scal_lmp, wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00557 
00558 
00559     double *nodes = (double*)calloc(P, sizeof(double));
00560     double *weights = (double*)calloc(P, sizeof(double));
00561     assert(nodes != NULL);
00562     assert(weights != NULL);
00563     flag_spherlaguerre_sampling(nodes, weights, R, P);
00564 
00565     flag_core_synthesis(f, flmp, nodes, P, L, P);
00566 
00567     free(nodes);
00568     free(weights);
00569     free(wav_lmp);
00570     free(scal_lmp);
00571     free(f_wav_lmp);
00572     free(f_scal_lmp);
00573     free(flmp);
00574 }
00575 
00595 void flaglet_axisym_wav_analysis_multires_real(double *f_wav, double *f_scal, const double *f, double R, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00596 {
00597     complex double *flmp;
00598     flag_core_allocate_flmn(&flmp, L, P);
00599     flag_core_analysis_real(flmp, f, R, L, P);
00600 
00601     double *wav_lmp, *scal_lmp;
00602     complex double *f_wav_lmp, *f_scal_lmp;
00603 
00604     flaglet_axisym_allocate_wav_lmp(&wav_lmp, &scal_lmp, B_l, B_p, L, P);
00605     flaglet_axisym_wav_lmp(wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00606 
00607     flaglet_axisym_allocate_f_wav_multires_lmp(&f_wav_lmp, &f_scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00608     flaglet_axisym_wav_analysis_multires_lmp(f_wav_lmp, f_scal_lmp, flmp, wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00609 
00610     free(wav_lmp);
00611     free(scal_lmp);
00612 
00613     double *nodes = (double*)calloc(P, sizeof(double));
00614     double *weights = (double*)calloc(P, sizeof(double));
00615     assert(nodes != NULL);
00616     assert(weights != NULL);
00617     flag_spherlaguerre_sampling(nodes, weights, R, P);
00618 
00619     int bandlimit_p, bandlimit_l, offset_lmp, offset, jl, jp;
00620     int J_l = s2let_j_max(L, B_l);
00621     int J_p = s2let_j_max(P, B_p);
00622     flag_core_synthesis_real(f_scal, f_scal_lmp, nodes, P, L, P);
00623 
00624     offset_lmp = 0;
00625     offset = 0;
00626     for (jp = J_min_p; jp <= J_p; jp++){
00627         bandlimit_p = MIN(s2let_bandlimit(B_p, jp), P);
00628         flag_spherlaguerre_sampling(nodes, weights, R, bandlimit_p);
00629         for (jl = J_min_l; jl <= J_l; jl++){
00630             bandlimit_l = MIN(s2let_bandlimit(B_l, jl), L);
00631             //offset_lmp = jp * (J_l + 1) * L * L * P   +  jl * L * L * P;
00632             //offset = jp * (J_l + 1) * L * (2*L-1) * (P)  +  jl * L * (2*L-1) * (P);
00633             flag_core_synthesis_real(f_wav + offset, f_wav_lmp + offset_lmp, nodes, bandlimit_p, bandlimit_l, bandlimit_p);
00634             offset_lmp += bandlimit_l * bandlimit_l * bandlimit_p ;
00635             offset += bandlimit_l * (2 * bandlimit_l - 1) * bandlimit_p;
00636         }
00637     }
00638 
00639     free(nodes);
00640     free(weights);
00641     free(f_wav_lmp);
00642     free(f_scal_lmp);
00643     free(flmp);
00644 }
00645 
00665 void flaglet_axisym_wav_synthesis_multires_real(double *f, const double *f_wav, const double *f_scal, double R, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00666 {
00667     complex double *f_wav_lmp, *f_scal_lmp;
00668     flaglet_axisym_allocate_f_wav_multires_lmp(&f_wav_lmp, &f_scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00669 
00670     int bandlimit_p, bandlimit_l, offset_lmp, offset, jl, jp;
00671     int J_l = s2let_j_max(L, B_l);
00672     int J_p = s2let_j_max(P, B_p);
00673     flag_core_analysis_real(f_scal_lmp, f_scal, R, L, P);
00674 
00675     offset_lmp = 0;
00676     offset = 0;
00677     for (jp = J_min_p; jp <= J_p; jp++){
00678         bandlimit_p = MIN(s2let_bandlimit(B_p, jp), P);
00679         for (jl = J_min_l; jl <= J_l; jl++){
00680             bandlimit_l = MIN(s2let_bandlimit(B_l, jl), L);
00681             //offset_lmp = jp * (J_l + 1) * L * L * P   +  jl * L * L * P;
00682             //offset = jp * (J_l + 1) * L * (2*L-1) * (P)   +  jl * L * (2*L-1) * (P);
00683             flag_core_analysis_real(f_wav_lmp + offset_lmp, f_wav + offset, R, bandlimit_l, bandlimit_p);
00684             offset_lmp += bandlimit_l * bandlimit_l * bandlimit_p ;
00685             offset += bandlimit_l * (2 * bandlimit_l - 1) * bandlimit_p;
00686         }
00687     }
00688 
00689     double *wav_lmp, *scal_lmp ;
00690     flaglet_axisym_allocate_wav_lmp(&wav_lmp, &scal_lmp, B_l, B_p, L, P);
00691     flaglet_axisym_wav_lmp(wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00692 
00693     complex double *flmp;
00694     flag_core_allocate_flmn(&flmp, L, P);
00695     flaglet_axisym_wav_synthesis_multires_lmp(flmp, f_wav_lmp, f_scal_lmp, wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00696     
00697     double *nodes = (double*)calloc(P, sizeof(double));
00698     double *weights = (double*)calloc(P, sizeof(double));
00699     assert(nodes != NULL);
00700     assert(weights != NULL);
00701     flag_spherlaguerre_sampling(nodes, weights, R, P);
00702 
00703     flag_core_synthesis_real(f, flmp, nodes, P, L, P);
00704 
00705     free(nodes);
00706     free(weights);
00707     free(wav_lmp);
00708     free(scal_lmp);
00709     free(f_wav_lmp);
00710     free(f_scal_lmp);
00711     free(flmp);
00712 }
00713 
00732 void flaglet_axisym_wav_analysis(complex double *f_wav, complex double *f_scal, const complex double *f, double R, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00733 {
00734     complex double *flmp;
00735     flag_core_allocate_flmn(&flmp, L, P);
00736     flag_core_analysis(flmp, f, R, L, P);
00737 
00738     double *wav_lmp, *scal_lmp;
00739     complex double *f_wav_lmp, *f_scal_lmp;
00740 
00741     flaglet_axisym_allocate_wav_lmp(&wav_lmp, &scal_lmp, B_l, B_p, L, P);
00742     flaglet_axisym_wav_lmp(wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00743 
00744     flaglet_axisym_allocate_f_wav_lmp(&f_wav_lmp, &f_scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00745     flaglet_axisym_wav_analysis_lmp(f_wav_lmp, f_scal_lmp, flmp, wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00746 
00747     free(wav_lmp);
00748     free(scal_lmp);
00749 
00750     double *nodes = (double*)calloc(P, sizeof(double));
00751     double *weights = (double*)calloc(P, sizeof(double));
00752     assert(nodes != NULL);
00753     assert(weights != NULL);
00754     flag_spherlaguerre_sampling(nodes, weights, R, P);
00755 
00756     int offset_lmp, offset, jl, jp;
00757     int J_l = s2let_j_max(L, B_l);
00758     int J_p = s2let_j_max(P, B_p);
00759     flag_core_synthesis(f_scal, f_scal_lmp, nodes, P, L, P);
00760 
00761     offset_lmp = 0;
00762     offset = 0;
00763     for (jp = J_min_p; jp <= J_p; jp++){
00764         for (jl = J_min_l; jl <= J_l; jl++){
00765             //offset_lmp = jp * (J_l + 1) * L * L * P   +  jl * L * L * P;
00766             //offset = jp * (J_l + 1) * L * (2*L-1) * (P)  +  jl * L * (2*L-1) * (P);
00767             flag_core_synthesis(f_wav + offset, f_wav_lmp + offset_lmp, nodes, P, L, P);
00768             offset_lmp += L * L * P ;
00769             offset += L * (2 * L - 1) * P;
00770         }
00771     }
00772 
00773     free(nodes);
00774     free(weights);
00775     free(f_wav_lmp);
00776     free(f_scal_lmp);
00777     free(flmp);
00778 }
00779 
00798 void flaglet_axisym_wav_synthesis(complex double *f, const complex double *f_wav, const complex double *f_scal, double R, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00799 {
00800     complex double *f_wav_lmp, *f_scal_lmp;
00801     flaglet_axisym_allocate_f_wav_lmp(&f_wav_lmp, &f_scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00802 
00803     int offset_lmp, offset, jl, jp;
00804     int J_l = s2let_j_max(L, B_l);
00805     int J_p = s2let_j_max(P, B_p);
00806     flag_core_analysis(f_scal_lmp, f_scal, R, L, P);
00807 
00808     offset_lmp = 0;
00809     offset = 0;
00810     for (jp = J_min_p; jp <= J_p; jp++){
00811         for (jl = J_min_l; jl <= J_l; jl++){
00812             //offset_lmp = jp * (J_l + 1) * L * L * P   +  jl * L * L * P;
00813             //offset = jp * (J_l + 1) * L * (2*L-1) * (P)   +  jl * L * (2*L-1) * (P);
00814             flag_core_analysis(f_wav_lmp + offset_lmp, f_wav + offset, R, L, P);
00815             offset_lmp += L * L * P ;
00816             offset += L * (2 * L - 1) * P;
00817         }
00818     }
00819 
00820     double *wav_lmp, *scal_lmp ;
00821     flaglet_axisym_allocate_wav_lmp(&wav_lmp, &scal_lmp, B_l, B_p, L, P);
00822     flaglet_axisym_wav_lmp(wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00823 
00824     complex double *flmp;
00825     flag_core_allocate_flmn(&flmp, L, P);
00826     flaglet_axisym_wav_synthesis_lmp(flmp, f_wav_lmp, f_scal_lmp, wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00827     
00828     double *nodes = (double*)calloc(P, sizeof(double));
00829     double *weights = (double*)calloc(P, sizeof(double));
00830     assert(nodes != NULL);
00831     assert(weights != NULL);
00832     flag_spherlaguerre_sampling(nodes, weights, R, P);
00833 
00834     flag_core_synthesis(f, flmp, nodes, P, L, P);
00835 
00836     free(nodes);
00837     free(weights);
00838     free(wav_lmp);
00839     free(scal_lmp);
00840     free(f_wav_lmp);
00841     free(f_scal_lmp);
00842     free(flmp);
00843 }
00844 
00864 void flaglet_axisym_wav_analysis_real(double *f_wav, double *f_scal, const double *f, double R, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00865 {
00866     complex double *flmp;
00867     flag_core_allocate_flmn(&flmp, L, P);
00868     flag_core_analysis_real(flmp, f, R, L, P);
00869 
00870     double *wav_lmp, *scal_lmp;
00871     complex double *f_wav_lmp, *f_scal_lmp;
00872 
00873     flaglet_axisym_allocate_wav_lmp(&wav_lmp, &scal_lmp, B_l, B_p, L, P);
00874     flaglet_axisym_wav_lmp(wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00875 
00876     flaglet_axisym_allocate_f_wav_lmp(&f_wav_lmp, &f_scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00877     flaglet_axisym_wav_analysis_lmp(f_wav_lmp, f_scal_lmp, flmp, wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00878 
00879     free(wav_lmp);
00880     free(scal_lmp);
00881 
00882     double *nodes = (double*)calloc(P, sizeof(double));
00883     double *weights = (double*)calloc(P, sizeof(double));
00884     assert(nodes != NULL);
00885     assert(weights != NULL);
00886     flag_spherlaguerre_sampling(nodes, weights, R, P);
00887 
00888     int offset_lmp, offset, jl, jp;
00889     int J_l = s2let_j_max(L, B_l);
00890     int J_p = s2let_j_max(P, B_p);
00891     flag_core_synthesis_real(f_scal, f_scal_lmp, nodes, P, L, P);
00892 
00893     offset_lmp = 0;
00894     offset = 0;
00895     for (jp = J_min_p; jp <= J_p; jp++){
00896         for (jl = J_min_l; jl <= J_l; jl++){
00897             //offset_lmp = jp * (J_l + 1) * L * L * P   +  jl * L * L * P;
00898             //offset = jp * (J_l + 1) * L * (2*L-1) * (P)  +  jl * L * (2*L-1) * (P);
00899             flag_core_synthesis_real(f_wav + offset, f_wav_lmp + offset_lmp, nodes, P, L, P);
00900             offset_lmp += L * L * P ;
00901             offset += L * (2 * L - 1) * P;
00902         }
00903     }
00904 
00905     free(nodes);
00906     free(weights);
00907     free(f_wav_lmp);
00908     free(f_scal_lmp);
00909     free(flmp);
00910 }
00911 
00931 void flaglet_axisym_wav_synthesis_real(double *f, const double *f_wav, const double *f_scal, double R, int B_l, int B_p, int L, int P, int J_min_l, int J_min_p)
00932 {
00933     complex double *f_wav_lmp, *f_scal_lmp;
00934     flaglet_axisym_allocate_f_wav_lmp(&f_wav_lmp, &f_scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00935 
00936     int offset_lmp, offset, jl, jp;
00937     int J_l = s2let_j_max(L, B_l);
00938     int J_p = s2let_j_max(P, B_p);
00939     flag_core_analysis_real(f_scal_lmp, f_scal, R, L, P);
00940 
00941     offset_lmp = 0;
00942     offset = 0;
00943     for (jp = J_min_p; jp <= J_p; jp++){
00944         for (jl = J_min_l; jl <= J_l; jl++){
00945             //offset_lmp = jp * (J_l + 1) * L * L * P   +  jl * L * L * P;
00946             //offset = jp * (J_l + 1) * L * (2*L-1) * (P)   +  jl * L * (2*L-1) * (P);
00947             flag_core_analysis_real(f_wav_lmp + offset_lmp, f_wav + offset, R, L, P);
00948             offset_lmp += L * L * P ;
00949             offset += L * (2 * L - 1) * P;
00950         }
00951     }
00952 
00953     double *wav_lmp, *scal_lmp ;
00954     flaglet_axisym_allocate_wav_lmp(&wav_lmp, &scal_lmp, B_l, B_p, L, P);
00955     flaglet_axisym_wav_lmp(wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00956 
00957     complex double *flmp;
00958     flag_core_allocate_flmn(&flmp, L, P);
00959     flaglet_axisym_wav_synthesis_lmp(flmp, f_wav_lmp, f_scal_lmp, wav_lmp, scal_lmp, B_l, B_p, L, P, J_min_l, J_min_p);
00960     
00961     double *nodes = (double*)calloc(P, sizeof(double));
00962     double *weights = (double*)calloc(P, sizeof(double));
00963     assert(nodes != NULL);
00964     assert(weights != NULL);
00965     flag_spherlaguerre_sampling(nodes, weights, R, P);
00966 
00967     flag_core_synthesis_real(f, flmp, nodes, P, L, P);
00968 
00969     free(nodes);
00970     free(weights);
00971     free(wav_lmp);
00972     free(scal_lmp);
00973     free(f_wav_lmp);
00974     free(f_scal_lmp);
00975     free(flmp);
00976 }
00977 
00992 int jjlmp2ind(int jl, int jp, int l, int m, int n, int J_l, int J_p, int L, int P)
00993 {
00994     return jp * (J_l + 1) * L * L * P   +  jl * L * L * P  +  n * L * L  +  l * l  + l +  m;
00995 }
00996 
01006 int lmp2ind(int l, int m, int n, int L)
01007 {
01008     return n * L * L  +  l * l  +  l  +  m ;
01009 }
01010 
01011