clc;
clear;
tic
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%首先设置系统静态参数%%%%%%%%
R=0.25; %线圈半径
N=30; %线圈匝数
I=1; %单匝电流
mu=4*pi*1e-7; %真空磁导率
SX=683.21; %X分量接收线圈等效面积
SY=717.34; %Y分量接收线圈等效面积
SZ=790.86; %Z分量接收线圈等效面积
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%设置目标体位置、倾角参数%%%%%%
x0=0.5; %设置目标体x坐标为测线中心
y0=0.5; %设置目标体y坐标测线中心
z0=-0.38; %设置目标体z坐标 -0.5m
fei=0*pi/4; %设置目标体fei角度
theta=1*pi/4; %设置目标体倾角参数
%%%%%目标体特征响应模型82mm迫击炮弹%%%%%
%%%沿主轴方向特征响应参数
k1=9.655e-4;
alfa1=25.6e-6;
beta1=0.8891;
gama1=5.242e-3;
%%%主轴垂直方向特征响应参数
k2=2.317e-5;
alfa2=-22.76e-6;
beta2=1.132;
gama2=2.616e-3;
%%%得到目标体的特性响应
t=5e-6:5e-6:20e-3; %衰减时间为20ms，采集时间间隔5us
L1_t=k1*(t+alfa1).^(-beta1).*exp(-t/gama1); %得到目标体主轴平行方向时域归一化响应
L2_t=k2*(t+alfa2).^(-beta2).*exp(-t/gama2); %得到目标体主轴垂直方向时域归一化响应
%%%特征响应抽道，抽道比例1.1406，抽35道
n=0:35;
c=1.1406.^n;
cc=20*c; %抽道起点为20
ccc=round(cc);
for i=1:35
L1_tt(i)=mean(L1_t(cc(i):cc(i+1)-1));
L2_tt(i)=mean(L2_t(cc(i):cc(i+1)-1));
end
L1_t=L1_tt;
L2_t=L2_tt;
%%%特征响应出图
% figure(1)
% loglog(1:35,L1_t,'*',1:35,L2_t,'s');grid on;
%计算特征响应矩阵B
for i=1:length(ccc)-1
B(:,:,i)=[L2_t(i),0,0;0,L2_t(i),0;0,0,L1_t(i)];
end
%欧拉旋转角
A=[cos(theta)*cos(fei),-sin(fei),sin(theta)*cos(fei);cos(theta)*sin(fei),cos(fei),sin(theta)*sin(fei);-sin(theta),0,cos(theta)];
%磁极化张量
for i=1:length(ccc)-1
M(:,:,i)=A*B(:,:,i)*A';
end
%磁极化方向向量
for i=1:length(ccc)-1
M1(:,:,i)=[M(1,1,i),M(1,2,i),M(1,3,i),M(2,2,i),M(2,3,i),M(3,3,i)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%设计二维计算条件%%%%%%%%%%%%
xL=1.0; %计算范围
dx=0.01; %x方向的分辨率
x=0:dx:xL; %x方向范围0到xL米
yL=1.0; %计算范围
dy=0.01; %y方向的分辨率
y=0:dy:yL; %y方向范围0到yL米
%%%%%%%%%%%%%%%%%计算当整个系统移动时，目标体处的一次场%%%%%%%%%%%%%%%%%
%%%系统坐标
[X,Y]=meshgrid(x,y); %生成系统的向量坐标r
r=zeros(length(y),length(x),3); %length(y)* length(x)个点，每个点的坐标包含x、y、z三个分量
r(:,:,1)=X;
r(:,:,2)=Y;
%%%目标体对发射系统的相对位置
R_r(:,:,1)=x0-r(:,:,1); %相对位置x分量
R_r(:,:,2)=y0-r(:,:,2); %相对位置y分量
R_r(:,:,3)=z0-r(:,:,3); %相对位置z分量
%%%计算每一个位置在目标体处产生的一次场
BP=zeros(length(y),length(x),3);
KK=mu*N*I/(2*pi); %系数
for i=1:length(y) %给定一个y坐标，算这一条测线上系统产生的一次场，再针对其他测线进行计算
xx=R_r(i,:,1); %相对位置的X分量
yy=y0-y(i); %相对位置的Y分量
%%%一次场计算过程见“核心文档——偶极子正演理论”式（21）
a=xx.^2+yy.^2+z0.^2+R^2;
b=xx*R;
c=yy*R;
d=z0*R;
M=1000; %将圆形线圈分为1000段
for n=0:M-1 %计算每一段圆弧在目标体处产生的一次场，总和即为线圈产生的一次场
fm=(a-2*b*cos(2*n*pi/M)-2*c*sin(2*n*pi/M)).^1.5;
BP(i,:,1)=BP(i,:,1)+d*cos(2*n*pi/M).*sin(pi/M)./fm; %一次场x分量
BP(i,:,2)=BP(i,:,2)+d*sin(2*n*pi/M).*sin(pi/M)./fm; %一次场y分量
BP(i,:,3)=BP(i,:,3)+(R.^2-c*sin(2*n*pi/M)-b*cos(2*n*pi/M)).*sin(pi/M)./fm; %一次场z分量
end
end
%%%一次场计算结果
BPX=BP(:,:,1)*KK;
BPY=BP(:,:,2)*KK;
BPZ=BP(:,:,3)*KK;
%%%一次场出图
% figure(2) %观察系统在不同位置时在目标体处产生的一次场
% subplot(1,3,1)
% contourf(X,Y,BPX);grid on; %一次场x分量
% subplot(1,3,2)
% contourf(X,Y,BPY);grid on; %一次场y分量
% subplot(1,3,3)
% contourf(X,Y,BPZ);grid on; %一次场z分量
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%二次场响应%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算目标体相对系统位置参量
r_R=-R_r;
RRX=r_R(:,:,1);
RRY=r_R(:,:,2);
RRZ=r_R(:,:,3);
%转换一次场的三个分量
HPX=BPX;
HPY=BPY;
HPZ=BPZ;
%定义参数矩阵中参量
for i=1:length(y)
for j=1:length(x)
RR(i,j)=sqrt(RRX(i,j).^2+RRY(i,j).^2+RRZ(i,j).^2);
k11(i,j)=(HPX(i,j)*(3*RRX(i,j)^2-(RR(i,j))^2))/(RR(i,j))^5;
k12(i,j)=(3*RRX(i,j)*(RRX(i,j)*HPY(i,j)+RRY(i,j)*HPX(i,j))-HPY(i,j)*(RR(i,j))*(RR(i,j)))/(RR(i,j))^5;
k13(i,j)=(3*RRX(i,j)*(RRX(i,j)*HPZ(i,j)+RRZ(i,j)*HPX(i,j))-HPZ(i,j)*(RR(i,j))*(RR(i,j)))/(RR(i,j))^5;
k14(i,j)=3*RRX(i,j)*RRY(i,j)*HPY(i,j)/(RR(i,j))^5;
k15(i,j)=(3*RRX(i,j)*(RRY(i,j)*HPZ(i,j)+RRZ(i,j)*HPY(i,j)))/(RR(i,j))^5;
k16(i,j)=3*RRX(i,j)*RRZ(i,j)*HPZ(i,j)/(RR(i,j))^5;
k21(i,j)=3*RRY(i,j)*RRX(i,j)*HPX(i,j)/(RR(i,j))^5;
k22(i,j)=(3*RRY(i,j)*(RRX(i,j)*HPY(i,j)+RRY(i,j)*HPX(i,j))-HPX(i,j)*(RR(i,j))*(RR(i,j)))/(RR(i,j))^5;
k23(i,j)=(3*RRY(i,j)*(RRX(i,j)*HPZ(i,j)+RRZ(i,j)*HPX(i,j)))/(RR(i,j))^5;
k24(i,j)=(HPY(i,j)*(3*RRY(i,j)^2-(RR(i,j))^2))/(RR(i,j))^5;
k25(i,j)=(3*RRY(i,j)*(RRY(i,j)*HPZ(i,j)+RRZ(i,j)*HPY(i,j))-HPZ(i,j)*(RR(i,j))*(RR(i,j)))/(RR(i,j))^5;
k26(i,j)=3*RRY(i,j)*RRZ(i,j)*HPZ(i,j)/(RR(i,j))^5;
k31(i,j)=3*RRZ(i,j)*RRX(i,j)*HPX(i,j)/(RR(i,j))^5;
k32(i,j)=(3*RRZ(i,j)*(RRX(i,j)*HPY(i,j)+RRY(i,j)*HPX(i,j)))/(RR(i,j))^5;
k33(i,j)=(3*RRZ(i,j)*(RRX(i,j)*HPZ(i,j)+RRZ(i,j)*HPX(i,j))-HPX(i,j)*(RR(i,j))*(RR(i,j)))/(RR(i,j))^5;
k34(i,j)=3*RRZ(i,j)*RRY(i,j)*HPY(i,j)/(RR(i,j))^5;
k35(i,j)=(3*RRZ(i,j)*(RRY(i,j)*HPZ(i,j)+RRZ(i,j)*HPY(i,j))-HPY(i,j)*(RR(i,j))*(RR(i,j)))/(RR(i,j))^5;
k36(i,j)=(HPZ(i,j)*(3*RRZ(i,j)^2-(RR(i,j))^2))/(RR(i,j))^5;
end
end
Q=1/(4*pi);
%计算二次电压V的三个分量
for i=1:length(ccc)-1
VX0(:,:,i)=Q*(k11*M1(1,1,i)+k12*M1(1,2,i)+k13*M1(1,3,i)+k14*M1(1,4,i)+k15*M1(1,5,i)+k16*M1(1,6,i));
end
for i=1:length(ccc)-1
VY0(:,:,i)=Q*(k21*M1(1,1,i)+k22*M1(1,2,i)+k23*M1(1,3,i)+k24*M1(1,4,i)+k25*M1(1,5,i)+k26*M1(1,6,i));
end
for i=1:length(ccc)-1
VZ0(:,:,i)=Q*(k31*M1(1,1,i)+k32*M1(1,2,i)+k33*M1(1,3,i)+k34*M1(1,4,i)+k35*M1(1,5,i)+k36*M1(1,6,i));
end
%计算水平分量响应
VXY0=sqrt(VX0.^2+VY0.^2);
%转换为mV
VX0=VX0*1e3;
VY0=VY0*1e3;
VXY0=VXY0*1e3;
VZ0=VZ0*1e3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%响应绘图分析%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%取第1道，观察早期响应
n=1;
figure(3)
subplot(2,4,1)
contourf(X,Y,VX0(:,:,n));title('X component');xlabel('x(m)');ylabel('y(m)');h=colorbar;set(get(h,'Title'),'string','mV');%grid on;
subplot(2,4,2)
contourf(X,Y,VY0(:,:,n));title('Y component');xlabel('x(m)');ylabel('y(m)');h=colorbar;set(get(h,'Title'),'string','mV');%grid on;
subplot(2,4,3)
contourf(X,Y,VXY0(:,:,n));title('Horizontal component');xlabel('x(m)');ylabel('y(m)');h=colorbar;set(get(h,'Title'),'string','mV');%grid on;
subplot(2,4,4)
contourf(X,Y,VZ0(:,:,n));title('Z component');xlabel('x(m)');ylabel('y(m)');h=colorbar;set(get(h,'Title'),'string','mV');%grid on;
n=35; %取第25道，观察晚期响应
subplot(2,4,5)
contourf(X,Y,VX0(:,:,n));title('X component');xlabel('x(m)');ylabel('y(m)');h=colorbar;set(get(h,'Title'),'string','mV');%grid on;
subplot(2,4,6)
contourf(X,Y,VY0(:,:,n));title('Y component');xlabel('x(m)');ylabel('y(m)');h=colorbar;set(get(h,'Title'),'string','mV');%grid on;
subplot(2,4,7)
contourf(X,Y,VXY0(:,:,n));title('Horizontal component');xlabel('x(m)');ylabel('y(m)');h=colorbar;set(get(h,'Title'),'string','mV');%grid on;
subplot(2,4,8)
contourf(X,Y,VZ0(:,:,n));title('Z component');xlabel('x(m)');ylabel('y(m)');h=colorbar;set(get(h,'Title'),'string','mV');
