hts.m

计算水声传播的快速场(FFP)程序

function p = HTS( kmin, deltak, stabil, r, G, option )

% Hankel transform of SCOOTER Green's function to produce pressure
% mbp, Dec. 2003

% G      = input k-space spectra, G( nrd, nt ) (destroyed by HTS)
% p      = output transmission loss, p( nrd, nr ).
% kmin   = minimum wavenumber
% deltak = spacing between wavenumbers
% stabil = amount integration has been moved off of the real axis.
% r      = vector of ranges for TL (must be evenly spaced)

% note: may be worth looking at Matlab tuning for fft (see fftw command)

if ndims( G ) > 2; error( 'Hankel transform called with an n-dimensional array, n>3' ); end;
if size( G, 2 ) == 1   % if G is a vector, expand it into a matrix with one row
    G = reshape( G, 1, length( G ) );
end

nk  = size( G, 2 );  % number of points in transform
nt2 = 2 * nk - 2;
k   = 0: deltak : ( nk - 1 ) * deltak;
r   = linspace( r( 1 ), r( end ), nt2 );
nrd = size( G, 1 );
ck  = stabil - i * kmin;
p   = zeros( nt2, nrd );

% scaling ...
if ( kmin > 0.5 * deltak )
   G( :, 1 ) =                       deltak   / sqrt( 2*pi ) * sqrt( kmin +             i*stabil ) * G( :, 1 );
else
   G( :, 1 ) =    0.5 * ( i*stabil + deltak ) / sqrt( 2*pi ) * sqrt(                    i*stabil ) * G( :, 1 );
end
G( :, 2:nk ) = scalecol( G( :, 2:nk ), deltak / sqrt( 2*pi ) * sqrt( kmin + k( 2:nk ) + i*stabil ) );
G( :, 1    ) = real( G( :, 1 ) );

switch option( 2: 2 )
    case 'P'   % positive spectrum
        Y = scalecol( G, exp( -i * r( 1 ) * k + i*pi/4 ) ).'; % Y transposed so that FFT done down columns
        Y = [ Y; zeros( nk-2, size( G, 1 ) ) ];
        Y = fft( Y, nt2 ).';    % exp (-IKX) transform; Y transposed so that each row is pressure vs. range
        p = scalecol( Y, exp ( ck * r ) );
    case 'N'   % negative spectrum
        Y = scalecol( G, exp( +i * r( 1 ) * k - i*pi/4 ) ).';
        Y = [ zeros( size( Y ) ); flipud( Y( 2:nk-1, : ) ) ];
        Y = fft( Y, nt2 ).';  % exp (+IKX) transform with normalization factor, nt
        p = scalecol( Y, exp ( -ck * r ) );
    case 'B'   % both positive and negative spectrum
        Y = scalecol( G, exp( -i * r( 1 ) * k + i*pi/4 ) ).'; % Y transposed so that FFT done down columns
        Y = [ Y; zeros( nk-2, size( G, 1 ) ) ];
        Y = fft( Y, nt2 ).';    % exp (-IKX) transform; Y transposed so that each row is pressure vs. range
        p = scalecol( Y, exp ( ck * r ) );
        
        Y = scalecol( G, exp( +i * r( 1 ) * k - i*pi/4 ) ).';
        Y = [ zeros( size( Y ) ); flipud( Y( 2:nk-1, : ) ) ];
        Y = fft( Y, nt2 ).';  % exp (+IKX) transform with normalization factor, nt
        p = p + scalecol( Y, exp ( -ck * r ) );    
end

% cylindrical spreading
if ( option(1:1) == 'R' )
    %ii = find( r < eps( max( abs( r ) ) ) ); % look for zeros (or near-zeros) in the range vector
    %r( ii ) = eps( max( abs( r ) ) );        % get rid of them
    r( r < eps( max( abs( r ) ) ) ) = eps( max( abs( r ) ) ); % get rid of zeros in the range vector
    p = scalecol( p, 1 ./ sqrt( abs( r ) ) );
end