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
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 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 % 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;