%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
clear
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 初始化

% 已知條件
wavelength=1.0; %  波長
a=wavelength;  % 令圓柱半徑等于波長
               % 為便于處理，令來波方向角為pi,入射角為0

% 常數(shù)定義
K=128; % 圓柱截面圓離散點數(shù)
dp=2.0*pi/K; % 離散角步長
del=2*pi*a/K; % 離散弧步長
k=2.0*pi/wavelength; % 波數(shù)
u0=4*pi*10^(-7);   % 真空磁導率
e0=1/(36*pi*10^9);   % 真空介電常數(shù)
ETA=(u0/e0)^(1/2);   % 自由空間波阻抗
gama=1.78107; % 歐拉常數(shù)
d=k*ETA/4.0;   %　系數(shù)
e=2.718;
% 坐標確定
for m=1:K
    pp(m)=(m-1)*dp+dp/2;
    x(m)=a*cos(pp(m));
    y(m)=a*sin(pp(m));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 主體

% 阻抗矩陣計算
for m=1:K      % 場點循環(huán)
    for n=1:K    % 源點循環(huán)
        if (n==m)   % 奇異點處理
            LMN(m,n)=k*ETA/4.0*del*(1-(2/pi)*i*(log(k*gama*del/(4.0*e))));   % 零階近似
%            
        else
            r=sqrt((x(m)-x(n))^2+(y(m)-y(n))^2);
            LMN(m,n)=k*ETA/4.0*del*besselh(0,2,k*r);   % 零階近似
%             
        end
    end
end  

% 電壓矩陣
for m=1:K;
    GM(m)=exp(-i*k*(a*cos(pp(m))+a*sin(pp(m))));
end

% % 表面電流分布矩陣計算
Jz=LMN\GM.';
disp(Jz)
figure(1)
plot(abs(Jz));
title('電流分布圖');