ita_sph_convSH.m 6.78 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204
function vector_out = ita_sph_convSH(vector_1,vector_2,MAX_degree)
%ITA_SPH_CONVSH - convolution of spherical harmonic coefficients
% function vector_out = ita_sph_convSH(vector1,vector2,MAX_degree)
%
% Function calculates the convolution of two sets of SH coefficients
% Input shoud be two vectors of size N^2 with the coefficients to be
% convoluted.
%
% MAX_degree allows for the user to define the up to which degree the
% convolution should be calculated. If not specified,
% MAX_degree = degree_1 + degree_2
%
% For velocity, when a given coefficient is calculated for the first time
% it is saved at a cell array inside the MAT-file "GauntValues",
% considerably increasing the calculation speed on the subsequent runs.
%
% The definition was taken from:
% [1] Driscoll and Healy, "Computing Fourier Transforms and Convolutions on the
% 2-D Sphere", Advances in Applied Mathematics, 15, 202-250, 1994
%
% [2] Messiah, A. "Quantum Mechanics", Courier Dover Publications, 1999
% in page 1057
%
% Bruno Masiero (bma@akustik.rwth-aachen.de)
% Institute of Technical Acoustics, RWTH Aachen, Germany,
% 19/08/2008

% <ITA-Toolbox>
% This file is part of the application SphericalHarmonics for the ITA-Toolbox. All rights reserved. 
% You can find the license for this m-file in the application folder. 
% </ITA-Toolbox>


[m,n] = size(vector_1);
if n > 1
    if m > 1
        error('Input should be a vector.')
    else
        vector_1 = vector_1(:);
    end
end
% assure that vector2 is a column vector
[m,n] = size(vector_2);
if n > 1
    if m > 1
        error('Input should be a vector.')
    else
        vector_2 = vector_2(:);
    end
end
% get degree and order of the inputs
size_vector_1 = length(vector_1); degree_1 = sqrt(size_vector_1)-1;
size_vector_2 = length(vector_2); degree_2 = sqrt(size_vector_2)-1;
if (rem(degree_1,1) ~= 0) || (rem(degree_2,1) ~= 0)
    error('Vectors should be of size N^2.');
end
% if vectors are not the same size, vector2 should be the largest one
if size_vector_1 > size_vector_2
    aux = vector_2;
    vector_2 = vector_1;
    vector_1 = aux;
    clear aux;
    
    size_vector_1 = length(vector_1); degree_1 = sqrt(size_vector_1)-1;
    size_vector_2 = length(vector_2); degree_2 = sqrt(size_vector_2)-1;
end
% generate degree and order vectors
[aux_degree_1,aux_order_1] = ita_sph_linear2degreeorder(1:size_vector_1);
[aux_degree_2,aux_order_2] = ita_sph_linear2degreeorder(1:size_vector_2);
% calculate the order and size of resulting vector
if nargin == 3
    degree_out = MAX_degree;
else
    degree_out = degree_1 + degree_2;
end
size_vector_out = (degree_out + 1)^2;
vector_out = zeros(size_vector_out,1);
%%%%%%%  Initiate Memory Matrix  %%%%%%%
if ~exist('GauntValues.mat')
    GauntCoefficients = cell(size_vector_1,size_vector_2);
else
    load('GauntValues.mat')
    [size_1,size_2] = size(GauntCoefficients);
    if  size_1 < size_vector_1
        GauntCoefficients = [GauntCoefficients;...
            cell(size_vector_1-size_1,size_2)];
        size_1 = size_vector_1;
    end
    if size_2 < size_vector_2
        GauntCoefficients = [GauntCoefficients ...
            cell(size_1,size_vector_2-size_2)];
    end
end
update = 0;
% Calculate Coefficients
for p = 1:size_vector_1
    for q = 1:size_vector_2
        l1 = aux_degree_1(p);
        l2 = aux_degree_2(q);
                      
        % Check if value was already calculated, remembering the diagonal
        % simetry of the cell array, since
        % (l1 m1 l2 m2 L M)(l1 0 l2 0 L 0) = (l2 m2 l1 m1 L M)(l2 0 l1 0 L 0).
        if q >= p
            G = GauntCoefficients{p,q};
        else
            G = GauntCoefficients{q,p};
        end
        
        % if not yet calculated, calculate it.
        if isempty(G)
            m1 = aux_order_1(p);
            m2 = aux_order_2(q);
            G = sparse(ita_sph_degreeorder2linear(l1+l2),1);
            
            % Y*Y = SUM(Y), see reference [2].
            for L = abs(l1-l2):l1+l2
                M = m1+m2;
                % if abs(M) > L, then Ylm is defined as 0.
                if abs(M) <= L
                    index = ita_sph_degreeorder2linear(L,M);
                    C1 = Wigner3j([l1 l2 L],[0 0 0]);
                    C2 = Wigner3j([l1 l2 L],[m1 m2 -M]);
                    cte = (-1)^M * sqrt((2*l1+1)*(2*l2+1)*(2*L+1)/(4*pi));
                    G(index) = cte*C1*C2;
                end
            end
            GauntCoefficients{p,q} = G;
            update = 1;
        end
        l12 = l1+l2;
        if l12 > degree_out
            l12 = degree_out;
        end
        index_l12 = ita_sph_degreeorder2linear(l12);
        
        vector_out(1:index_l12) = vector_out(1:index_l12) + ...
                (G(1:index_l12)*vector_1(p)*vector_2(q));
    end
end
if update
    save('GauntValues.mat', 'GauntCoefficients')
end
        
function wigner = Wigner3j( j123, m123 )
% Compute the Wigner 3j symbol using the Racah formula. 
%
% W = Wigner3j( J123, M123 ) 
%
% J123 = [J1, J2, J3], with the condition:
%        |Ji - Jj| <= Jk <= (Ji + Jj)    (i,j,k are permutations of 1,2,3)
% M123 = [M1, M2, M3], with the conditions:
%        |Mi| <= Ji    (i = 1,2,3)
%        M1 + M2 + M3 = 0
% All Ji and Mi have to be half integers (correspondingly).
% 
% Reference: 
% Wigner 3j-Symbol entry of Eric Weinstein's Mathworld:
% http://mathworld.wolfram.com/Wigner3j-Symbol.html
%
% Inspired by Wigner3j.m by David Terr, Raytheon, 6-17-04
%  (available at www.mathworks.com/matlabcentral/fileexchange).
%
% By Kobi Kraus, Technion, 25-6-08.
% Modified by BMA on 19/08/2008
j1 = j123(1); j2 = j123(2); j3 = j123(3);
m1 = m123(1); m2 = m123(2); m3 = m123(3);
% Input error checking
if any( j123 < 0 )
    wigner = 0; return%error( 'The j must be non-negative' )
elseif any( rem( [j123, m123], 0.5 ) )
    wigner = 0; return%error( 'All arguments must be integers or half-integers' );
elseif any( rem( (j123 - m123), 1 ) | ( abs( m123 ) > j123 ) )
    wigner = 0; return%error( 'j123 and m123 do not match' );
elseif ( j3 > (j1 + j2) ) || ( j3 < abs(j1 - j2) )
    wigner = 0; return%error( 'j3 is out of bounds' );
elseif m1 + m2 + m3 ~= 0
    wigner = 0; return%error( 'm3 does not match m1 + m2' );
end
% Simple common case
if ~any( m123 ) && rem( sum( j123 ), 2 ) % m1 = m2 = m3 = 0 & j1 + j2 + j3 is odd
    wigner = 0;
    return
end
% Calculation
t1 = j2 - m1 - j3;
t2 = j1 + m2 - j3;
t3 = j1 + j2 - j3;
t4 = j1 - m1;
t5 = j2 + m2;
tmin = max( 0,  max( t1, t2 ) );
tmax = min( t3, min( t4, t5 ) );
t = tmin : tmax;
wigner = sum( (-1).^t .* exp( -ones(1,6) *...
    gammaln( [t; t-t1; t-t2; t3-t; t4-t; t5-t] +1 ) + ...
    gammaln( [j1+j2+j3+1, j1+j2-j3, j1-j2+j3, -j1+j2+j3, j1+m1, j1-m1, j2+m2, j2-m2, j3+m3, j3-m3] +1 )...
     * [-1; ones(9,1)] * 0.5 ) ) * (-1)^( j1-j2-m3 );
         
% Warnings
% if isnan( wigner )
%     warning( 'Wigner3J is NaN!' )
% elseif isinf( wigner )
%     warning( 'Wigner3J is Inf!' )
% end