24小时热门版块排行榜

返回列表

当前主题已经存档。

【有奖交流】积极回复本帖子，参与交流，就有机会分得作者 lzjuxmch 的 10 个金币

lzjuxmch

铁虫 (小有名气)

应助: 0 (幼儿园)
金币: 299.3
帖子: 112
在线: 44.3小时
虫号: 490252
注册: 2008-01-03
专业: 凝聚态物性 II ：电子结构

[交流] 【求助】hankel变换翻译成C语言

本人急需hankel变换的C程序，目前只能找到matlab中hankel变换的原代码，但是本人不懂matlab，请各位帮个忙能否帮把下列代码翻译成C或者C++语言？提供下列程序的思路也可。万分感谢！

%[H,I]=ht(h,r,k,I)
%-----------------
%
%Hankel transform of order 0.
%
%Input:
% h    Signal h(r)
% r    Radial positions [m]       {0:numel(h)-1}
% k    Spatial frequencies [rad/m] {pi/numel(h)*(0:numel(h)-1)}
% I    Integration kernel ?       {default}
%
%Output:
% H    Spectrum H(k)
% I    Integration kernel
%
%
% ?  If the integration kernel is missing, it is
%    recomputed from the Bessel functions (slow).
%

%    Marcel Leutenegger ?June 2006
%
function [H,I]=ht(h,r,k,I)
if sum(size(h) > 1) > 1
error('Signal must be a vector.');
end
if nargin < 2 | isempty(r)
r=0:numel(h)-1;
else
[r,w]=sort(r(

.');
h=h(w);
end
if nargin < 3 | isempty(k)
k=pi/numel(h)*(0:numel(h)-1);
end
if nargin < 4 | isempty(I)
r=[(r(2:end) + r(1:end-1))/2 r(end)];
I=2*pi./k(

*r.*besselj(1,k(

*r);
I(k == 0,

=pi*r.*r;
I=I - [zeros(numel(k),1) I(:,1:end-1)];
elseif exist('w','var')
I=I(:,w);
end
H=reshape(I*h(

,size(k));

[ Last edited by mingdong on 2009-4-16 at 17:53 ]

回复此楼