util.c

Go to the documentation of this file.

 
#include "s2kit/util.h"
 
#include <math.h>
#include <stdlib.h>
 
void inline ComplexMult(const double x, const double y, const double u, const double v,
                        double* real_result, double* imag_result) {
    *real_result = x * u - y * v;
    *imag_result = x * v - y * u;
}
 
int IndexOfHarmonicCoeff(const int m, const int l, const int bw) {
    int bigL = bw - 1;
 
    if (m >= 0)
        return (m * (bigL + 1) - ((m * (m - 1)) / 2) + (l - m));
 
    return (((bigL * (bigL + 3)) / 2) + 1 + ((bigL + m) * (bigL + m + 1) / 2) + (l - abs(m)));
}
 
void TransMult(double* rdatacoeffs, double* idatacoeffs, double* rfiltercoeffs, double* ifiltercoeffs, double* rres,
               double* ires, const int bw) {
    double* rdptr = rdatacoeffs;
    double* idptr = idatacoeffs;
    double* rrptr = rres;
    double* irptr = ires;
 
    for (int m = 0; m < bw; ++m) {
        for (int l = m; l < bw; ++l) {
            ComplexMult(rfiltercoeffs[l], ifiltercoeffs[l], rdptr[l - m], idptr[l - m], rrptr + l - m, irptr + l - m);
 
            rrptr[l - m] *= sqrt(4. * M_PI / (2. * l + 1.));
            irptr[l - m] *= sqrt(4. * M_PI / (2. * l + 1.));
        }
        rdptr += bw - m;
        idptr += bw - m;
        rrptr += bw - m;
        irptr += bw - m;
    }
 
    int size = 2 * bw;
    for (int m = bw + 1; m < size; ++m) {
        for (int l = size - m; l < bw; ++l) {
            ComplexMult(rfiltercoeffs[l], ifiltercoeffs[l], rdptr[l - size + m], idptr[l - size + m],
                        rrptr + l - size + m, irptr + l - size + m);
 
            rrptr[l - size + m] *= sqrt(4. * M_PI / (2. * l + 1.));
            irptr[l - size + m] *= sqrt(4. * M_PI / (2. * l + 1.));
        }
 
        rdptr += m - bw;
        idptr += m - bw;
        rrptr += m - bw;
        irptr += m - bw;
    }
}