I need to do a spherical harmonic analysis
if you have a routine that produces
expansion coefficients Cnm Snm
please help
thx,
A
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
>
>
>
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: