Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

Spherical Harmonic Cnm Snm

126 views
Skip to first unread message

Alex S

unread,
Jan 19, 2003, 8:24:32 PM1/19/03
to
hello,

I need to do a spherical harmonic analysis
if you have a routine that produces
expansion coefficients Cnm Snm

please help
thx,
A


Alex S

unread,
Jan 20, 2003, 12:37:33 PM1/20/03
to
Thanks, Timo
I am a non-mathematician and have problems
to implement integration in Matlab
(not mentioning the other two methods)
How do I obtain Ynm in Matlab?
and then integrate over sphere?
I am trying to do a Real ( Cnm and Snm)
spherical harmonic expansion,
I used LEGENDRE function to obtain Pnm like this

Pn = sqrt(2*n(i)+1)*legendre(n(i),sin(lats),'sch');
Pnm = reshape(p(m(i)+1,:,:),nlat,nlon);

and then i am trying to integrate by the following...
-- is this a right way to do it (please see below the complete
listing)?----

> Snm(i) = (1/(length(n)*4*pi)) *
sum(sum(DATA.*sin((m(i))*lons.*Pnm));
....
> Pnm = Pnm.*cos( m(i)*lons);
....
> Cnm(i) = (DATA * Pnm')./(Pnm * Pnm');
>
When calculating real spherical harmonics I need to track somehow
even/odd
m,n?

I'll greatly appreciate your help

thanks,
A
>
>
> "Timo Nieminen" <ti...@physics.uq.edu.au> wrote in message
> news:Pine.LNX.4.33L2.03012...@localhost.localdomain...
> > On Sun, 19 Jan 2003, Alex S wrote:
> >
> > > I am trying to figure out a spherical harmonic routine to derive
> > > Cnm, Snm coefficients,
> >
> > You have (at least) three choices:
> >
> > (1) Use an integral transform:
> >
> > Cnm = int( data .* conj(Ynm) )/int( Ynm .* conj(Ynm) );
> >
> > where the integrals are over the sphere. Note that the denominator = 1,
> > since the spherical harmonics are orthonormal. Note that if you are
using
> > real spherical harmonics, the complex conjugate is unnecessary.
> >
> > (2) Use point-matching:
> >
> > Decide to what maximum n you want to find the expansion to. Then write
> > down a system of linear equations relating the expansion coefficients to
> > the data. Solve this system. Works best if your linear system is
> > overdetermined.
> >
> > (3) Use a fast spherical transform. See eg Lesur, Geophys. J. Int. 139,
> > 547-555 (1999)
> >
> > --
> > Timo Nieminen - Home page: http://www.physics.uq.edu.au/people/nieminen/
> > Shrine to Spirits:
http://www.users.bigpond.com/timo_nieminen/spirits.html
> >
> hi,
> I am trying to figure out a spherical harmonic routine to derive
> Cnm, Snm coefficients,
>
> please help
> thx,
> A
>
>
>
> function [c,s,m,n] = SHT(DATA,m,n)
>
> hgt = DATA;
> res = 4; nlon = 360/res+1;nlat = 180/res+1;
> lons = linspace(0,360,nlon)/180*pi;
> lats = linspace(90, -90,nlat)/180*pi;
> [lons,lats] = meshgrid(lons,lats);
>
> hgt1 = reshape(hgt,[1 prod(size(hgt))]);
>
> for i = 1:length(n)
>
> if (m(i) == 0)
> p = sqrt(2*n(i)+1)*legendre(n(i),sin(lats),'sch');
> end
> if (n(i) == 0)
> pp = p;
> else
> pp = reshape(p(m(i)+1,:,:),nlat,nlon);
>
> end
> s(i) = (1/(length(n)*4*pi)) * sum(sum(hgt.*sin((m(i))*lons.*pp));
> pp = pp.*cos((m(i))*lons);
> pp = reshape(pp, [1 prod(size(pp))]);
> c(i) = (hgt1 * pp')./(pp * pp');
> end
>
>
>


Christian Gruber

unread,
Jan 23, 2003, 7:09:57 AM1/23/03
to
Hi,
A very simple manner is to just summarize over
the discrete function values of a 1 degree grid of phi, lam.
The functions, here a potential comes e.g. on a 181x360 grid.
In each step the entire latitude range gets processed.
Note that the coefficients A,B for cnm, snm are both stored in a
diagonal scheme in one matrix.


function [Cnm] = decomp2D(V, degree);
% field = decomp2D(V, degree)
% V - the discrete potentials
%
% CG, TU-Berlin 01/02

constants;

[m,n] = size(V);
step = 360/max(m,n);
if mod(360,step) ~= 0
error('wrong stepsize');
end


rho = 180/pi;
Th = [0:step:180];

Cnm = zeros(degree+1,degree+1);
pind = Th*step+1;


for L = 0:step:359
for n = 0:degree
for m = 0:n

lind = L*step+1;

% Alm
Cnm(n+1,m+1) = Cnm(n+1,m+1) + ...
sum(V(pind, lind) .* plm(n,m,Th)' .* sin(Th/rho)') *
cos(m*L/rho) * (step/rho)^2;
if L == 359
Cnm(n+1,m+1) = Cnm(n+1,m+1) * 1/4/pi;
end
if n == 0
continue;
end
if m == 0
continue;
end
% Blm
Cnm(m,n+1) = Cnm(m,n+1) + ...
sum(V(pind, lind) .* plm(n,m,Th)' .* sin(Th/rho)') *
sin(m*L/rho) * (step/rho)^2;
if L == 359
Cnm(m,n+1) = Cnm(m,n+1) * 1/4/pi;
end
end
end
end


Alex S schrieb:

0 new messages