liov.m

时域矩量法解线天线辐射问题和散射场问题

%%%%完纯导体地架空线雷电感应过电压计算%%%%%
clear all
tic


%定义通用变量
mu=4*pi*10^(-7);
c=3*10^8;
epsilon0=8.85e-12;



% l=1;
% a=0.002;
% h=0.5;


l=10;
a=0.05;
h=30;

rr=100;



M=31;
N=1e2;
dx=l/M;
dt=dx/c;
center=(M+1)/2;
%预设结果存储矩阵
I=zeros(M,N);
F=zeros(M,N);
E=zeros(M,N);



K=zeros(M,M);
K1=zeros(M,M);
K2=zeros(M,M);


PA=zeros(M,N);
PA1=zeros(M,N);
PA2=zeros(M,N);


A=zeros(M,N);
A1=zeros(M,N);
A2=zeros(M,N);

%%%%%%%%%%%%%%%%%%%%%%%计算中间系数%%%%%%%%%%%%%%%%%%%%%
for m=1:M
    for k=1:M
        za=(m-k)*dx+dx/2;
        zb=(m-k)*dx-dx/2;
        K1(m,k)=mu/4/pi*(log(za+(za^2+a^2)^0.5)-log(zb+(zb^2+a^2)^0.5));
        K2(m,k)=mu/4/pi*(log(za+(za^2+4*h^2)^0.5)-log(zb+(zb^2+4*h^2)^0.5));
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%激励为雷电辐射场%%%%%%%%%%%%
%%%从雷电通道到架空线做垂线，设垂足为cz点，设cz=center%%%

% cz=center;
% for m=1:M
%     r=(rr^2+(abs(m-cz)*dx)^2)^0.5;
%     for n=1:N
%        tend=n*dt;
%        re=Er_time(r,tend);
%        E(m,n)=re;
%      end
% end

%%%%%%%%激励为高斯脉冲函数%%%%%%%%%%%%%%%%%%%%%
t0=4e-9;
tao=4e-8;
for m=1:M
    for n=1:N
        E(m,n)=exp(-4*pi*(n*dt-t0)^2/tao^2);
    end
end
% figure(1)
% hold on
% for m=1:M
% plot(E(m,:))
% end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%激励调试区%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%雷电辐射场%%%%%%%%%%%%%%%%%%%%%%%%%%

% dett=3e-7;
% te=30e-6;
% tt=0e-6:dett:te;
% E1=zeros(1,length(tt));
%
% for n=1:length(tt)
%       tend=tt(n);
% %        re=Er_time_xps;
%        re=Er_time;
%        E1(n)=re;
% end


%%%%%%%%高斯脉冲函数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% t0=2e-9;
% tao=4e-9;
% % dtt=3e-10;
% n=1:1:N;
t=n*dt;
% E1=exp(-4*pi*((t-t0)/tao).^2);
plot(t,E(1,:));

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



%%%%%%各点入射场强偏导数值%%%%%%%%%%%%%%%%%%%%%%
for m=1:M
    for n=1:(N-1)
        F(m,n)=(E(m,n+1)-E(m,n))/((c^2)*(dt));
    end
    F(m,N)=F(m,N-1);
end
n=3;
while n<=N
    for m=1:M
        for k=1:M
            tt1=n-abs(m-k);
            if tt1>0
                A1(m,n)=A1(m,n)+I(k,tt1)*K1(m,k);
                if k~=m
                    PA1(m,n)=PA1(m,n)+I(k,tt1)*K1(m,k);
                end
            end
            tt2=fix((n*dt)-(((m-k)*dx)^2+4*h^2)^0.5/c);
            if tt2>0
                A2(m,n)=A2(m,n)+I(k,tt2)*K2(m,k);
                if k~=m
                    PA2(m,n)=PA2(m,n)+I(k,tt2)*K1(m,k);
                end
            end
        end
        A(m,n)=A1(m,n)+A2(m,n);
        PA(m,n)=PA1(m,n)+PA2(m,n);
    end
    for m=2:M-1
        I(m,n)=1/K1(m,m)*(-PA(m,n)+2*A(m,n-1)-A(m,n-2)-(c*dt)^2*F(m,n-1)...
            +(c*dt/dx)^2*(A(m+1,n-1)-2*A(m,n-1)+A(m-1,n-1)));
    end
    n=n+1;
end

for n=1:N
    x(n)=n*dt;
    y(n)=I(center,n);
end
plot(x,y),grid

toc