# v_sphrharm

## PURPOSE

V_SPHRHARM forward and inverse spherical harmonic transform

## SYNOPSIS

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

## DESCRIPTION

    V_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);
                    v_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 (includes neither pole) [default]
                          for invertibility, ne>=2N+1 for order N
                   'U' uniform inclination grid:  e=(0:ne-1)*pi/ne (includes North pole only)
                          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=v_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=v_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]=v_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:

- v_cblabel V_CBLABEL add a label to a colorbar C=(L,H)

This function is called by:

## SUBFUNCTIONS

- function [ueo,weo,lgpo]=sphrharp(gr,ne,nsp)

## SOURCE CODE

    0001 function [u,v,w]=v_sphrharm(m,a,b,c,d)
    0002 %V_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 %                v_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 (includes neither pole) [default]
    0034 %                      for invertibility, ne>=2N+1 for order N
    0035 %               'U' uniform inclination grid:  e=(0:ne-1)*pi/ne (includes North pole only)
    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=v_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=v_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]=v_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=v_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-v_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(v_sphrharm('fur',v_sphrharm('iur',h,10,10),no)-h))
    0132 % same but complex: no=4; h=rand(1,(no+1)^2); max(abs(v_sphrharm('fu',v_sphrharm('iu',h,10,10),no)-h))
    0133 % same but gaussian grid: no=4; h=rand(1,(no+1)^2); max(abs(v_sphrharm('fg',v_sphrharm('ig',h,6,10),no)-h))
    0134
    0135 %      Copyright (C) Mike Brookes 2009
    0136 %      Version: $Id: v_sphrharm.m 10865 2018-09-21 17:22:45Z 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 Lesser General Public License as published by
    0144 %   the Free Software Foundation; either version 3 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 Lesser General Public License for more details.
    0151 %
    0152 %   You can obtain a copy of the GNU Lesser General Public License from
    0153 %   https://www.gnu.org/licenses/ .
    0154 %
    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             %             v_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                 v_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; % vector of elevations
    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;