Need help: Spherical harmonic integration

3 次查看（过去 30 天）
Yousuf Khan
Yousuf Khan 2021-12-15
编辑： Yousuf Khan 2022-1-3
Let's say I have a combination of spherical harmonics e.g.
Ytotal=aY_4^2+bY_6^2+cY_8^2 ([1]
So, because spherical harmonics are an orthogonal basis, I can say:
eq2) b=∫Y_6^2 * Y_total ( [2]
Now, I have a function that gives me a spherical harmonic (The famous spharm4)
function varargout=plm2spec(lmcosi,norma,in,ot)
% [sdl,l,bta,lfit,logy,logpm]=PLM2SPEC(lmcosi,norma,in,ot)
%
% Calculates the power spectrum of real spherical harmonic
% sine and cosine coefficients contained in the matrix 'lmcosi' which is
% of the standard form that can be plotted by PLM2XYZ and PLOTPLM.
%
% INPUT:
%
% lmcosi Spherical harmonic coefficients [l m Ccos Csin]
% norma 1 multiplication by (l+1)
% This gives the mean-square value of the
% gradient of a potential in Schmidt-harmonics
% 2 division by (2*l+1) [default]
% This gives the proper power spectral density
% as we've come to know it
% 3 none, i.e. a scaling factor of 1
% in Index to minimum degree to consider for the spectral fit [defaulted]
% ot Index to maximum degree to consider for the spectral fit [defaulted]
%
% OUTPUT:
%
% sdl Spectral density: energy per degree
% l Degree
% bta Spectral slope of loglog(l,sdl)
% lfit,logy Spectral line plot given by loglog(lfit,logy)
% logpm Error on spectral line plot given by
% loglog(lfit,logpm)
%
% SEE ALSO: MTVAR
%
% EXAMPLE:
%
% [sdl,l,bta,lfit,logy,logpm]=plm2spec(fralmanac('EGM96'));
%
% SEE ALSO: ACTSPEC
%
% The normalization by (2l+1) is what's required when the spherical
% harmonics are normalized to 4pi. See DT p. 858. A "delta"-function then
% retains a flat spectrum. See Dahlen and Simons 2008.
% See papers by Hipkin 2001, Kaula 1967, Lowes 1966, 1974, Nagata 1965
% (Lowes, JGR 71(8), 2179 [1966])
% (Nagata, JGeomagGeoel 17, 153-155 [1965])
%
% Last modified by fjsimons-at-alum.mit.edu, 03/18/2020
defval('norma',2)
lmin=lmcosi(1);
lmax=lmcosi(end,1);
pin=0;
for l=lmin:lmax
clm=shcos(lmcosi,l);
slm=shsin(lmcosi,l);
pin=pin+1;
sdl(pin)=clm(:)'*clm(:)+slm(:)'*slm(:);
end
switch norma
case 1
normfac=(lmin:lmax)+1;
case 2
normfac=1./(2*(lmin:lmax)+1);
case 3
normfac=1;
disp('Not further normalized')
otherwise
error('No valid normalization specified')
end
% disp(sprintf('Normalization %i',norma))
sdl=normfac.*sdl;
sdl=sdl(:);
% Figure out the range over which to make the fit
l=lmin:lmax;
l=l(:);
if lmin==0
defval('in',3);
elseif lmin==1
defval('in',2);
else
defval('in',1);
end
defval('ot',lmax)
lfit=l(in:ot);
if nargout>=3
% Calculate spectral slope
[bt,E]=polyfit(log10(lfit),log10(sdl(in:ot)),1);
bta=bt(1);
[logy,loge]=polyval(bt,log10(lfit),E);
logy=10.^logy;
logpm=[logy./(10.^loge) logy.*(10.^loge)];
else
[bta,lfit,logy,logpm]=deal(NaN);
end
which gives a spherical harmonic matrix
First, I want to check if the Y_6^2 is normalized (the integral should be equal to zero) using trapz. How can a 3D integral be done with trapz function using spherical coordinates?
After that, I want to find the b coefficient using eq(2), but I still can’t understand how to use the trapz function correctly and multiply the matrixs.
TL;DR how to use trapz to calculate a triple integral on a spherical harmonic?
Any help would be appreciated.

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by