%  MFIE.m
%  本程序用于磁場積分方程求解圓柱散射
%  編程人    沙威(Wei Sha) 安徽大學(Anhui University) ws108@ahu.edu.cn

clear;clc;
%%%% 1.常數

f=300*10^8;                 %  TM 波頻率
v=3*10^8;                   %  波速
anglemin=pi/128;            %  最小離散角度步長
amplitude=1;                %  電場振幅

wavelength=v/f;             %  真空波長
k0=2*pi/wavelength;         %  真空波數
a=4*wavelength;             %  圓柱半徑
count=2*pi/anglemin;        %  選配點數
circlelength=a*anglemin;    %  弧長

%%%%  2.坐標定位

%  圓柱中心建立坐標系
%  圓弧坐標

anglevalue=-anglemin/2;  %  角度初值

for m=1:count;
    anglevalue=anglevalue+anglemin;  %  角度遞增
    point(m)=anglevalue;             %  角度賦值
end;



%%%%  3.矩陣賦值

for m=1:count;                     %  場循環
    for n=1:count;                 %  源循環
        
        field=point(m);            %  場點位置
        source=point(n);           %  源點位置
        
        if (m==n)
            L(m,n)=-1/2;           %  奇異點處理
        else
            quad_value=2*k0*a*abs( sin( (field-source)/2 ) );
            L(m,n)=(j*k0/4)*circlelength*(besselj(1,quad_value)-j*bessely(1,quad_value))*...
                (abs( sin( (field-source)/2 ) ));
        end;
    end;
end;

for t=1:count;  %  電壓矩陣
    y(t)=-amplitude*exp(-j*a*cos(point(t))*k0);
end;



%%%%  4.求逆

result=L\y.';
figure(1);
plot(abs(result));
title('矩量法電流');



%%%%  5.rcs計算

t=0;
for s=0:pi/180:pi
    t=t+1;
    sum=0;
    for r=1:count;
        nr=cos(point(r)-s);  %  N.R,即遠場與源點方向
        sum=sum+nr*result(r,1)*circlelength*...
            exp(j*k0*a*cos(point(r))*cos(s)+j*k0*a*sin(point(r))*sin(s));
    end;
    rcs1(t)=(k0/4)*abs(sum).^2;
end;

rcs=10.*log10(rcs1/wavelength);  %  波長歸一化

figure(2);
plot(rcs);
title('雙站雷達散射截面');



% %%%%  6.解析解電流
% 
% for m=1:count;
%     sum=0;
%     for n=-30:30;
%         t=(1/2)*(besselh(n+1,2,k0*a)-besselh(n-1,2,k0*a));  %  bessel函數求導
%         sum=sum+((-j*2*amplitude)/(pi*k0*a))*j^(-n)*exp(j*n*(point(m)))/t;
%     end;
%     current(m)=sum;
% end;
% 
% figure(3);
% title('電流對比');
% plot(abs(current),'r');
% title('電流分布對比');
% hold on;
% plot(abs(result),'bx');
% legend('解析','矩量')
% 
% 
% 
% %%%%  7.解析解RCS
% 
% RR=10.^6;  %  足夠大的R
% t=0;
% for angle=0:pi/60:pi;
%     sum=0;
%     t=t+1;
%     for n=-30:30
%         cons2=(1/2)*(besselh(n+1,2,k0*a)-besselh(n-1,2,k0*a));
%         cons1=(1/2)*(besselj(n+1,k0*a)-besselj(n-1,k0*a));
%         A=cons1/cons2;
%         sum=sum+amplitude*j.^(-n)*A*besselh(n,2,k0*RR)*exp(j*n*angle);
%     end;
%     es(t)=sum;
% end;
% 
% RCS_ANA=2*pi*RR*abs(es).^2;
% RCS=10*log10(RCS_ANA/wavelength);
 
% figure(4);
% title('雷達散射截面對比');
% plot(RCS,'r');
% hold on;
% plot(rcs,'bx');
% legend('解析解','矩量解');