Home > voicebox > sphrharm.m

sphrharm

PURPOSE ^

SPHRHARM forward and inverse spherical harmonic transform

SYNOPSIS ^

function [u,v,w]=sphrharm(m,a,b,c,d)

DESCRIPTION ^

SPHRHARM  forward and inverse spherical harmonic transform

 Usage: (1) y=('f',n,x)      % Calculate complex transform of spatial data x up to order n

        (2) y=('fr',n,x)     % Calculate real transform of spatial data x(ne,na) up to order n
                             % x is specified on a uniform grid with inclination e=(0.5:ne-0.5)*pi/ne
                             % and azimuth a=(0:na-1)*2*pi/na. The North pole has inclination 0 and
                             % azimuths increase going East.

        (3) y=('fd',n,e,a,v) % Calculate transform of impulse at (e(i),a(i)) of height v(i) up to order n
                             % e(i), a(i) and v(i) should be the same dimension; v defaults to 1.

        (4) x=('i',y,ne,na)  % Calculate spatial data from spherical harmonics y on
                             % a uniform grid with ne inclinations and na azimuths

        (5) [e,a,w]=('cg',ne,na)   % Calculate the inclinations, azimuths and
                                   % quadrature weights of a Gaussian sampling grid

        (6) n=2;             % illustrate real basis functions upto order 2
            for i=0:n
              for j=-i:i
                subplot(n+1,2*n+1,i*(2*n+1)+j+n+1);
                sphrharm('irp',[zeros(1,i^2+i+j) 1],25,25);
              end
            end

  Inputs:  m string specifying various options:
               'f' forward transform
               'i' inverse transform
               'c' calculate the coordinates of the sampling grid [default if f or i omitted]
               'r' real spherical harmonics instead of complex
               'u' uniform inclination grid:  e=(0.5:ne-0.5)*pi/ne [default]
                      for invertibility, ne>=2N+1 for order N
               'U' uniform inclination grid:  e=(0:ne-1)*pi/ne
                      for invertibility, ne>=2N+1 for order N
               'g' gaussian inclination grid (non-uniform but fewer samples needed)
                      for invertibility, ne>=N+1 for order N
               'a' arbitrary (specified by user - inverse transform only)
               'd' delta function [forward transform only]
               'p' plot result
               'P' plot with colourbar

           The remaining inputs depend on the mode specified:
               'f'  a         order of transform
                    b(ne,na)  input data array on the chosen inclination-azimuth sampling grid
               'fd' a         order of transform
                    b(k)      inclinations of k delta functions
                    c(k)      azimuths of k delta functions
                    d(k)      amplitudes of k delta functions [default=1]
               'i'  a()       spherical harmonics as a single row vector in the order:
                                 (0,0),(1,-1),(1,0),(1,1),(2,-2),(2,-1),...
                              To access as a 2-D harmonic: Y(n,m)=a(n^2+n+m+1) where n=0:N, m=-n:n
                    b         number of inclination values in output grid
                    c         number of azimuth values in output grid
               'c'  a         number of inclination values in output grid
                    b         number of azimuth values in output grid

 Outputs:  u output data according to transform direction:
                'f': spherical harmonics as a single row vector in the order:
                        (0,0),(1,-1),(1,0),(1,1),(2,-2),(2,-1),...
                     To access as a 2-D harmonic: Y(n,m)=u(n^2+n+m+1) where n=0:N, m=-n:n
                'i': u(ne,na) is spatial data sampled on an azimuth-inclination grid
                     with ne inclination points (in 0:pi) and na azimuth points (in 0:2*pi)
                'c': u(ne) gives inclination grid positions with 0 being the North pole
           v(na) gives azimuth grid positions
           w(ne) gives the quadrature weights used in the transform

 Suppose f(e,a) is a complex-valued function defined on the surface of the
 sphere (0<=e<=pi, 0<=a<2pi) where e=inclination=pi/2-elevation and
 a=azimuth. (0,*) is the North pole, (pi/2,0) is on the equator near
 Ghana and (pi/2,pi/2) is on the equator near India.

 We can approximate f(e,a) using complex spherical harmonics as
 f(e,a) = sum(F(n,m)*u(n,m,e,a)) where the sum
 is over n=0:N,m=-n:n giving a total of (N+1)^2 coefficients F(n,m).

 If f(e,a) happens to be real-valued, then we can transform the
 coefficients using G(n,0)=F(n,0) and, for m>0,
 [G(n,m); G(n,-m)]=sqrt(0.5)*[1 1; 1i -1i]*[F(n,m); F(n,-m)]
 to give the real spherical harmonic coefficients G(n,m).

 The basis functions u(n,m,e,a) are orthonormal under the inner product
 <u(e,a),v(e,a)> = integral(u(e,a)*v'(e,a)*sin(e)*de*da).

  Minimum spatial grid for invertibility:
     Gaussian grid: ne >= N+1,  na >= 2N+1
     Uniform grid:  ne >= 2N+1, na >= 2N+1
     An inverse transform followed by a forward transform will restore the
     original coefficient values provided that the sampling grid is large
     enough. The advantage of the Gaussian grid is that it can be smaller.

   Data formats:
     (1) Spatial Data: x=sphrharm('i',y,ne,na)
            x(1:ne,1:na) is spatial data sampled on an azimuth-inclination grid
                         with ne inclination points (in 0:pi) and
                         na azimuth points (in 0:2*pi)
     (2) Spherical harmonics: y=sphrharm('f',n,x)
            y(1:(n+1)^2)  spherical harmonics as a single row vector
                          in the order: (0,0),(1,-1),(1,0),(1,1),(2,-2),(2,-1),...
                          To access as a 2-D harmonic, use
                          Y(n,m)=y(n^2+n+m+1)   where n=0:N, m=-i:i
     (3) Sampling Grid:  [e,a,w]=sphrharm('c',ne,na)
            e(1:ne)   monotonically increasing inclination values (0<=e<=pi) where 0 is the North pole
            a(1:na)   monotonically increasing azimuth values (0<=a<2pi)
            w(1:ne)   Quadrature weights

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [u,v,w]=sphrharm(m,a,b,c,d)
0002 %SPHRHARM  forward and inverse spherical harmonic transform
0003 %
0004 % Usage: (1) y=('f',n,x)      % Calculate complex transform of spatial data x up to order n
0005 %
0006 %        (2) y=('fr',n,x)     % Calculate real transform of spatial data x(ne,na) up to order n
0007 %                             % x is specified on a uniform grid with inclination e=(0.5:ne-0.5)*pi/ne
0008 %                             % and azimuth a=(0:na-1)*2*pi/na. The North pole has inclination 0 and
0009 %                             % azimuths increase going East.
0010 %
0011 %        (3) y=('fd',n,e,a,v) % Calculate transform of impulse at (e(i),a(i)) of height v(i) up to order n
0012 %                             % e(i), a(i) and v(i) should be the same dimension; v defaults to 1.
0013 %
0014 %        (4) x=('i',y,ne,na)  % Calculate spatial data from spherical harmonics y on
0015 %                             % a uniform grid with ne inclinations and na azimuths
0016 %
0017 %        (5) [e,a,w]=('cg',ne,na)   % Calculate the inclinations, azimuths and
0018 %                                   % quadrature weights of a Gaussian sampling grid
0019 %
0020 %        (6) n=2;             % illustrate real basis functions upto order 2
0021 %            for i=0:n
0022 %              for j=-i:i
0023 %                subplot(n+1,2*n+1,i*(2*n+1)+j+n+1);
0024 %                sphrharm('irp',[zeros(1,i^2+i+j) 1],25,25);
0025 %              end
0026 %            end
0027 %
0028 %  Inputs:  m string specifying various options:
0029 %               'f' forward transform
0030 %               'i' inverse transform
0031 %               'c' calculate the coordinates of the sampling grid [default if f or i omitted]
0032 %               'r' real spherical harmonics instead of complex
0033 %               'u' uniform inclination grid:  e=(0.5:ne-0.5)*pi/ne [default]
0034 %                      for invertibility, ne>=2N+1 for order N
0035 %               'U' uniform inclination grid:  e=(0:ne-1)*pi/ne
0036 %                      for invertibility, ne>=2N+1 for order N
0037 %               'g' gaussian inclination grid (non-uniform but fewer samples needed)
0038 %                      for invertibility, ne>=N+1 for order N
0039 %               'a' arbitrary (specified by user - inverse transform only)
0040 %               'd' delta function [forward transform only]
0041 %               'p' plot result
0042 %               'P' plot with colourbar
0043 %
0044 %           The remaining inputs depend on the mode specified:
0045 %               'f'  a         order of transform
0046 %                    b(ne,na)  input data array on the chosen inclination-azimuth sampling grid
0047 %               'fd' a         order of transform
0048 %                    b(k)      inclinations of k delta functions
0049 %                    c(k)      azimuths of k delta functions
0050 %                    d(k)      amplitudes of k delta functions [default=1]
0051 %               'i'  a()       spherical harmonics as a single row vector in the order:
0052 %                                 (0,0),(1,-1),(1,0),(1,1),(2,-2),(2,-1),...
0053 %                              To access as a 2-D harmonic: Y(n,m)=a(n^2+n+m+1) where n=0:N, m=-n:n
0054 %                    b         number of inclination values in output grid
0055 %                    c         number of azimuth values in output grid
0056 %               'c'  a         number of inclination values in output grid
0057 %                    b         number of azimuth values in output grid
0058 %
0059 % Outputs:  u output data according to transform direction:
0060 %                'f': spherical harmonics as a single row vector in the order:
0061 %                        (0,0),(1,-1),(1,0),(1,1),(2,-2),(2,-1),...
0062 %                     To access as a 2-D harmonic: Y(n,m)=u(n^2+n+m+1) where n=0:N, m=-n:n
0063 %                'i': u(ne,na) is spatial data sampled on an azimuth-inclination grid
0064 %                     with ne inclination points (in 0:pi) and na azimuth points (in 0:2*pi)
0065 %                'c': u(ne) gives inclination grid positions with 0 being the North pole
0066 %           v(na) gives azimuth grid positions
0067 %           w(ne) gives the quadrature weights used in the transform
0068 %
0069 % Suppose f(e,a) is a complex-valued function defined on the surface of the
0070 % sphere (0<=e<=pi, 0<=a<2pi) where e=inclination=pi/2-elevation and
0071 % a=azimuth. (0,*) is the North pole, (pi/2,0) is on the equator near
0072 % Ghana and (pi/2,pi/2) is on the equator near India.
0073 %
0074 % We can approximate f(e,a) using complex spherical harmonics as
0075 % f(e,a) = sum(F(n,m)*u(n,m,e,a)) where the sum
0076 % is over n=0:N,m=-n:n giving a total of (N+1)^2 coefficients F(n,m).
0077 %
0078 % If f(e,a) happens to be real-valued, then we can transform the
0079 % coefficients using G(n,0)=F(n,0) and, for m>0,
0080 % [G(n,m); G(n,-m)]=sqrt(0.5)*[1 1; 1i -1i]*[F(n,m); F(n,-m)]
0081 % to give the real spherical harmonic coefficients G(n,m).
0082 %
0083 % The basis functions u(n,m,e,a) are orthonormal under the inner product
0084 % <u(e,a),v(e,a)> = integral(u(e,a)*v'(e,a)*sin(e)*de*da).
0085 %
0086 %  Minimum spatial grid for invertibility:
0087 %     Gaussian grid: ne >= N+1,  na >= 2N+1
0088 %     Uniform grid:  ne >= 2N+1, na >= 2N+1
0089 %     An inverse transform followed by a forward transform will restore the
0090 %     original coefficient values provided that the sampling grid is large
0091 %     enough. The advantage of the Gaussian grid is that it can be smaller.
0092 %
0093 %   Data formats:
0094 %     (1) Spatial Data: x=sphrharm('i',y,ne,na)
0095 %            x(1:ne,1:na) is spatial data sampled on an azimuth-inclination grid
0096 %                         with ne inclination points (in 0:pi) and
0097 %                         na azimuth points (in 0:2*pi)
0098 %     (2) Spherical harmonics: y=sphrharm('f',n,x)
0099 %            y(1:(n+1)^2)  spherical harmonics as a single row vector
0100 %                          in the order: (0,0),(1,-1),(1,0),(1,1),(2,-2),(2,-1),...
0101 %                          To access as a 2-D harmonic, use
0102 %                          Y(n,m)=y(n^2+n+m+1)   where n=0:N, m=-i:i
0103 %     (3) Sampling Grid:  [e,a,w]=sphrharm('c',ne,na)
0104 %            e(1:ne)   monotonically increasing inclination values (0<=e<=pi) where 0 is the North pole
0105 %            a(1:na)   monotonically increasing azimuth values (0<=a<2pi)
0106 %            w(1:ne)   Quadrature weights
0107 %
0108 
0109 % future options
0110 %    direction: [m=rotation matrix]
0111 %    transform: [h=complex harmonics but only the positive ones specified]
0112 %    grid       d=delta function [forward transform only]
0113 %               [s=sparse azimuth]
0114 %               [z=include both north and south pole]
0115 %
0116 % bugs:
0117 %   (1) we could save space and time by taking advantage of symmetry of legendre polynomials
0118 %   (2) the number of points in the inclination (elevation) grid must, for now, be even if the 'U' option is chosen
0119 %   (3) should ensure that the negative nyquist azimuth frequency is not used
0120 %   (4) save time by manipulating only the necessary 2*m columns of the da matrix
0121 %   (5) should make non-existant coefficients black in plots
0122 %   (6) using 'surf' for plots adds an offset in azimuth and elevation
0123 %       because colours correspond to vertices not faces
0124 %   (7) the normalization for mode 'fd' seems incorrect
0125 %   (8) mode 'fd' should allow multiple impulses to be summed
0126 
0127 % errors:
0128 
0129 % tests:
0130 % check for inverse transform n=4; m=4; [ve,va]=spvals(8,8,n,m); ve*va-sphrharm('iur',[zeros(1,n^2+n+m) 1],8,8)
0131 % check inverse followed by forward: no=4; h=rand(1,(no+1)^2); max(abs(sphrharm('fur',sphrharm('iur',h,10,10),no)-h))
0132 % same but complex: no=4; h=rand(1,(no+1)^2); max(abs(sphrharm('fu',sphrharm('iu',h,10,10),no)-h))
0133 % same but gaussian grid: no=4; h=rand(1,(no+1)^2); max(abs(sphrharm('fg',sphrharm('ig',h,6,10),no)-h))
0134 
0135 %      Copyright (C) Mike Brookes 2009
0136 %      Version: $Id: sphrharm.m 2467 2012-11-02 08:12:41Z dmb $
0137 %
0138 %   VOICEBOX is a MATLAB toolbox for speech processing.
0139 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
0140 %
0141 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0142 %   This program is free software; you can redistribute it and/or modify
0143 %   it under the terms of the GNU General Public License as published by
0144 %   the Free Software Foundation; either version 2 of the License, or
0145 %   (at your option) any later version.
0146 %
0147 %   This program is distributed in the hope that it will be useful,
0148 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0149 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0150 %   GNU General Public License for more details.
0151 %
0152 %   You can obtain a copy of the GNU General Public License from
0153 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0154 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
0155 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0156 
0157 % decode option string
0158 mv='c u ';       % mv(1)=transform, mv(2)=real/complex, mv(3)=grid, mv(4)=plot
0159 if ~nargout
0160     mv(4)='P';
0161 end
0162 mc='ficruUgadpP';
0163 mi=[1 1 1 2 3 3 3 3 3 4 4];
0164 for i=1:length(m)
0165     j=find(m(i)==mc);
0166     if ~isempty(j)
0167         mv(mi(j))=m(i);
0168     end
0169 end
0170 
0171 switch mv(1)
0172     case 'f' % forward transform [sp]=('f',order,data)
0173         if mv(3)~='d'
0174             % input data has elevations down the columns and azimuths along the rows
0175             % output data is in the form:
0176             [ne,na]=size(b);
0177             da=fft(b,[],2)*sqrt(pi)/na;
0178             if mv(2)=='r' % only actually need to do this for min(a,floor(na/2)) values (but nyquist is tricky)
0179                 ix=2:ceil(na/2);
0180                 iy=na:-1:na+2-ix(end);
0181                 da(:,ix)=da(:,ix)+da(:,iy);
0182                 da(:,iy)=1i*(da(:,ix)-2*da(:,iy));  % note
0183                 da(:,1)=da(:,1)*sqrt(2);
0184             else
0185                 da=da*sqrt(2);
0186             end
0187             [ue,we,lgp]=sphrharp(mv(3),ne,a);
0188             da=da.*repmat(we',1,na);
0189             u=zeros(1,(a+1)^2);
0190             i=0;
0191             for nn=0:a
0192                 % we could vectorize this to avoid the inner loop
0193                 % we should ensure the the negative nyquist value is not used
0194                 for mm=-nn:nn
0195                     i=i+1;
0196                     u(i)=lgp(1+nn*(nn+1)/2+abs(mm),:)*da(:,1+mod(mm,na));
0197                 end
0198             end
0199         else       % forward transform of impulses [sp]=('fd',order,e,a,value)
0200             if nargin<5
0201                 d=ones(size(b));
0202             end
0203             [ue,we,lgp]=sphrharp('a',b(:),a);
0204             uu=zeros(1,(a+1)^2);     % reserve space for spherical harmonic coefficients
0205             u=uu;
0206             for j=1:numel(b)
0207                 i=0;
0208                 if mv(2)=='r'
0209                     for nn=0:a
0210                         % we could vectorize this to avoid the inner loop
0211                         % we should ensure the the negative nyquist value is not used
0212                         i=i+1;
0213                         uu(i+nn)=lgp(1+nn*(nn+1)/2,1);
0214                         for mm=1:nn
0215                             uu(i+nn-mm)=lgp(1+nn*(nn+1)/2+abs(mm),j)*sin(mm*c(j))*sqrt(2);
0216                             uu(i+nn+mm)=lgp(1+nn*(nn+1)/2+abs(mm),j)*cos(mm*c(j))*sqrt(2);
0217                         end
0218                         i=i+2*nn;
0219                     end
0220                 else
0221                     for nn=0:a
0222                         % we could vectorize this to avoid the inner loop
0223                         % we should ensure the the negative nyquist value is not used
0224                         for mm=-nn:nn
0225                             i=i+1;
0226                             uu(i)=lgp(1+nn*(nn+1)/2+abs(mm),j)*exp(-1i*mm*c(j));
0227                         end
0228                     end
0229                 end
0230                 u=u+d(j)*uu/sqrt(2*pi);
0231             end
0232         end
0233     case 'i' % inverse transform [data]=('i',sp,nincl,nazim)
0234         %or [data]=('a',sp,e,a) where e is a list of elevations and a is either a corresponding list of azimuths
0235         % or else a cell array contining several azimuths for each inclination
0236         % input data is sp=[(0,0) (1,-1) (1,0) (1,1) (2,-2) (2,-1) (2,0) ... ]
0237         % length is (n+1)^2 where n is the order
0238         % output data is an array [ne,na]
0239         nsp=ceil(sqrt(length(a)))-1;
0240         [ue,we,lgp]=sphrharp(mv(3),b,nsp);
0241         if mv(3)=='a'
0242             na=nsp*2+1;
0243             ne=length(b);
0244         else
0245             na=c;
0246             ne=b;
0247         end
0248         i=0;
0249         da=zeros(ne,na);
0250         for nn=0:nsp
0251             % this could be vectorized somewhat to speed it up
0252             for mm=-nn:nn
0253                 i=i+1;
0254                 if i>length(a)
0255                     break;
0256                 end
0257                 da(:,1+mod(mm,na))=da(:,1+mod(mm,na))+a(i)*lgp(1+nn*(nn+1)/2+abs(mm),:)';
0258             end
0259         end
0260         if na>1 && mv(2)=='r' % convert real to complex only actually need to do this for min(b,floor(na/2)) values (but nyquist is a bit tricky)
0261             ix=2:ceil(na/2);
0262             iy=na:-1:na+2-ix(end);
0263             da(:,iy)=(da(:,ix)+1i*da(:,iy))*0.5;
0264             da(:,ix)=da(:,ix)-da(:,iy);
0265             da(:,1)=da(:,1)/sqrt(2);
0266         else
0267             da=da/sqrt(2);
0268         end
0269         if mv(3)=='a' % do a slow fft
0270             if iscell(c)
0271                 u{ne,1}=[]; % reserve space for output cell array
0272                 for i=1:ne
0273                     ai=c{i};
0274                     ui=repmat(da(i,1),1,length(ai));
0275                     for j=1:nsp
0276                         exj=exp(1i*j*ai(:).'); % we could vectorize this more
0277                         ui=ui+da(i,j+1).'.*exj+da(i,na+1-j).'.*conj(exj);
0278                     end
0279                     u{i}=ui/sqrt(pi);
0280                 end
0281             else
0282                 u=da(:,1);
0283                 for j=1:nsp
0284                     exj=exp(1i*j*c(:)); % we could vectorize this more
0285                     u=u+da(:,j+1).*exj+da(:,na+1-j).*conj(exj);
0286                 end
0287                 u=u/sqrt(pi);
0288                 if mv(2)=='r' && isreal(a)
0289                     u=real(u);
0290                 end
0291             end
0292         else
0293             u=ifft(da,[],2)*na/sqrt(pi); % could put the scale factor 1/sqrt(pi) earlier
0294             if mv(2)=='r' && isreal(a)
0295                 u=real(u);
0296             end
0297         end
0298     case 'c' % just output coordinates [inclination,azim,weights]=('c',nincl,nazim)
0299         % if m!='g', then order, a, must be even
0300         [u,w]=sphrharp(mv(3),a,0);
0301         if nargin<3
0302             b=ceil(a/1+(mv(3)=='g'));
0303         end
0304         v=(0:b-1)*2*pi/b;
0305 end
0306 if mv(4)~=' '
0307     switch mv(1)
0308         case 'f'
0309             if mv(2)=='r'
0310                 ua=u;
0311                 tit='Real Coefficients';
0312             else
0313                 ua=abs(u);
0314                 tit='Complex Coefficient Magnitudes';
0315             end
0316             nu=length(ua);
0317             no=ceil(sqrt(nu))-1;
0318             yi=zeros(no,2*no+1);
0319             for i=0:no
0320                 for j=-i:i
0321                     iy=i^2+i+j+1;
0322                     if iy<=nu
0323                         yi(i+1,j+no+1)=ua(iy);
0324                     end
0325                 end
0326             end
0327             imagesc(-no:no,0:no,yi);
0328             axis 'xy';
0329             if mv(4)=='P'
0330                 colorbar;
0331             end
0332             xlabel('Harmonic');
0333             ylabel('Order');
0334             title(tit);
0335         case 'i'            % [data]=('i',sp,nincl,nazim)
0336             vv=(0:na)*2*pi/na;  % azimuth array
0337             gr=sin(ue)';
0338             surf(gr*cos(vv),gr*sin(vv),repmat(cos(ue)',1,na+1),real(u(:,[1:na 1])));
0339             axis equal;
0340             xlabel('X');
0341             ylabel('Y');
0342             zlabel('Z');
0343             if mv(4)=='P'
0344                 colorbar;
0345             end
0346             %             cblabel('Legendre Weight');
0347             %             title(sprintf('Sampling Grid: %s, %d, %d',mv(3),length(u),na-1));
0348         case 'c'            % [inclination,azim,weights]=('c',nincl,nazim)
0349             vv=[v v(1)];    % replicate initial azimuth point
0350             na=length(vv);
0351             gr=sin(u)';
0352             mesh(gr*cos(vv),gr*sin(vv),repmat(cos(u)',1,na),repmat(w',1,na));
0353             axis equal;
0354             xlabel('X');
0355             ylabel('Y');
0356             zlabel('Z');
0357             if mv(4)=='P'
0358                 colorbar;
0359                 cblabel('Quadrature Weight');
0360             end
0361             title(sprintf('Sampling Grid: %s, %d, %d',mv(3),length(u),na-1));
0362     end
0363 end
0364 
0365 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0366 
0367 function [ueo,weo,lgpo]=sphrharp(gr,ne,nsp)
0368 % calculate persistent variables for grid points and legendre polynomials
0369 % we recalculate if transform order or inclination grid size or type are changed
0370 % gr = grid type:
0371 %       'u'=uniform inclination starting at pi/2n (default)
0372 %       'U'=uniform but starting at north pole
0373 %       'g'=gaussian inclination
0374 %       'a'=arbitrary (specified by user - inverse transform only)
0375 % ne = number of inclination values
0376 % nsp = maximum order
0377 % Outputs:
0378 %    ueo(ne)     vector containing the inclination values
0379 %    weo(ne)     vector containing the quadrature weights
0380 %    lgpo((nsp+1)*(nsp+2)/2,ne)  evaluates the Schmidt seminormalized
0381 %                                associated Legendre function at each value
0382 %                                of cos(ue) for n=0:nsp, m=0:n.
0383 persistent gr0 lgp ne0 ue we nsp0
0384 if isempty(ne0)
0385     ne0=-1;
0386     nsp0=-1;
0387     gr0='a';
0388 end
0389 if gr=='a'
0390     ue=ne;
0391     ne=length(ue);
0392     ne0=-1;
0393     gr0=gr;
0394     nsp0=-1; % delete previous legendre polynomials
0395     lgp=[];
0396     we=[];
0397 elseif gr~=gr0 || ne~=ne0
0398     if gr=='g'
0399         r = 1:ne-1;
0400         r = r ./ sqrt(4*r.^2 - 1);
0401         p = zeros( ne, ne );
0402         p( 2 : ne+1 : ne*(ne-1) ) = r;
0403         p( ne+1 : ne+1 : ne^2-1 ) = r;
0404         [q, ue] = eig(p);
0405         [ue, k] = sort(diag(ue));
0406         ue=acos(-ue)';
0407         we = 2*q(1,k).^2;
0408     elseif gr=='U'
0409         if rem(ne,2)
0410             error('inclination grid size must be even when using ''U'' option');
0411         end
0412         ue=(0:ne-1)*pi/ne;
0413         xx=zeros(1,ne);
0414         ah=ne/2;
0415         xx(1:ah)=(1:2:ne).^(-1);
0416         we=-4*sin(ue).*imag(fft(xx).*exp(-1i*ue))/ne;
0417     else % default is m='u'
0418         ue=(1:2:2*ne)*pi*0.5/ne;
0419         vq=(ne-abs(ne+1-2*(1:ne))).^(-1).*exp(-1i*(ue+0.5*pi));
0420         we=(-2*sin(ue).*real(fft(vq).*exp(-1i*(0:ne-1)*pi/ne))/ne);
0421     end
0422     gr0=gr;
0423     ne0=ne;
0424     nsp0=-1; % delete previous legendre polynomials
0425     lgp=[];
0426 end
0427 if nsp>nsp0
0428     lgp((nsp+1)*(nsp+2)/2,ne)=0; % reserve space
0429     for i=nsp0+1:nsp
0430         lgp(1+i*(i+1)/2:(i+1)*(i+2)/2,:)=legendre(i,cos(ue),'sch')*sqrt(0.5*i+0.25);
0431         lgp(1+i*(i+1)/2,:)=lgp(1+i*(i+1)/2,:)*sqrt(2);
0432     end
0433     nsp0=nsp;
0434 end
0435 lgpo=lgp;
0436 weo=we;
0437 ueo=ue;

Generated on Tue 10-Oct-2017 08:30:10 by m2html © 2003