function LAYMT2TE
%层状均匀介质的二维正演程序，整体法剖分，TE模式
%TE模式
% MT2TE
KONGQILAY=5;%空气测层数
NXTE=60;
NYTE=60;
NPSUMTE=(NXTE+1)*(NYTE+1);%节点总数
NETE=NXTE*NYTE;%单元总数
xTE=zeros(NPSUMTE,1);
a=100;
h=100;
LLTE(1:NXTE)=a*ones(NXTE,1)';%横向网格长度，等距
BBTE(1:NYTE)=h*ones(NYTE,1)';%纵向网格宽度，等距,要考虑到去肤深度的因素，即503*sqre(rho/f)
ELBTE(1:NETE,1)=1:NETE;
%f=[0.001 0.0025 0.005 0.0075 0.0125 0.025 0.05 0.1 0.125 0.25 0.5 1 2 4 8 16 32 64 128 256 512 1024 2048];
%f=[0.001 0.005 0.01 0.02 0.05 0.1 0.125 0.2 0.25 0.5 1 1.25 2 2.5 5 10 12.5 20 25 50 100 500 1000];
f=[0.03125 0.0625 0.125 0.25 0.5 1 2 4 8 16 32 64 128 256 512 1024 2048];
PaTE=[100000000,100,500,1500];
for i=1:NXTE %单元编号
for j=1:NYTE
m=(i-1)*NYTE+j;
ELBTE(m,2)=LLTE(i);%单元的横向长度
ELBTE(m,3)=BBTE(j);%单元的纵向宽度
end
end
%%%%%%%%%%%%%%%%%%%%%%%%% 存放节点的X,Y编号‘MTJIEDIANBH’ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for IX=1:NXTE
for IY=1:NYTE
N=(IX-1)*NYTE+IY;%单元编号
N1=N+(IX-1);%每个单元的第一个节点编号
BHXY(1,N)=N1;%每个单元的第一个节点编号
BHXY(2,N)=N1+1;%每个单元的第二个节点编号
BHXY(3,N)=N1+NYTE+2;%每个单元的第三个节点编号
BHXY(4,N)=N1++NYTE+1;%每个单元的第四个节点编号
end
end
%%%%%%%%%%%%%%%%%%%%%%%%% 模型参数 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% for i=1:NXTE %空气层电阻率
% for j=1:KONGQILAY
% m=(i-1)*NYTE+j;
% rhoTE(m)=PaTE(1);
% end
% end
% for i=1:NXTE %第一层电阻率
% for j=KONGQILAY+1:NYTE
% m=(i-1)*NYTE+j;
% rhoTE(m)=PaTE(2);
% end
% end
% for i=21:40 %第二层电阻率
% for j=16:35
% m=(i-1)*NYTE+j;
% rhoTE(m)=PaTE(3);
% end
% end
for i=1:NXTE %空气层电阻率
for j=1:KONGQILAY
m=(i-1)*NYTE+j;
rhoTE(m)=PaTE(1);
end
end
for i=1:NXTE %1
for j=KONGQILAY+1:NYTE
m=(i-1)*NYTE+j;
rhoTE(m)=PaTE(2);%100
end
end
for i=15:29 %异常体
for j=9:23
m=(i-1)*NYTE+j;
rhoTE(m)=PaTE(3);%1000
end
end
for i=32:46 %异常体
for j=11:25
m=(i-1)*NYTE+j;
rhoTE(m)=PaTE(4);%500
end
end
for ff=1:size(f,2)
K1TE=zeros(NPSUMTE); %???????
K2TE=zeros(NPSUMTE); %???????
K3TE=zeros(NPSUMTE); %???????
PTE=zeros(NPSUMTE,1); %????????
u=(4e-7)*pi;%电、磁性参数
w=2*pi*f(ff);
TAOTE=1/(sqrt(-1)*w*u);
LMTE=1./rhoTE;
%%%%%%%%%%%%%%%%%%%% 形成K1e 及其扩展矩阵K1TE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for m=1:NETE
A=a/(h*6);
B=h/(a*6);
K1e=B*[2 1 -1 -2;1 2 -2 -1;-1 -2 2 1;-2 -1 1 2]+A*[2 -2 -1 1;-2 2 1 -1;-1 1 2 -2;1 -1 -2 2];
for j=1:4 %给每个单元的节点乘上K的系数
NJ=BHXY(j,m);
for k=1:4
NK=BHXY(k,m);
K1TE(NJ,NK)=K1TE(NJ,NK)+K1e(j,k)*TAOTE;%
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%% 形成K2e及其扩展矩阵K2TE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for m=1:NETE
K2e=[4 2 1 2;2 4 2 1;1 2 4 2;2 1 2 4];
for j=1:4 %给每个单元的节点乘上K的系数
NJ=BHXY(j,m);
for k=1:4
NK=BHXY(k,m);
K2TE(NJ,NK)=K2TE(NJ,NK)+K2e(j,k)*LMTE(m)*ELBTE(m,2)*ELBTE(m,3)/36;%??????????????? 如何处理重合叠加的Kij值
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 形成K3e及其扩展矩阵K3TE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for h1=NYTE:NYTE:NETE
i=BHXY(1,h1); j=BHXY(2,h1);
kn=sqrt(-sqrt(-1)*w*u./rhoTE(h1));
mk=TAOTE*kn*ELBTE(h1,2)/6;
Kjj=2*mk;
K3TE(j,j)=K3TE(j,j)+Kjj;
K3TE(i,i)=K3TE(i,i)+Kjj;
Kji=1*mk;
K3TE(j,i)=K3TE(j,i)+Kji;
K3TE(i,j)=K3TE(i,j)+Kji;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 组装总体刚度矩阵 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
vTE=sparse(K1TE-K2TE+K3TE);
%size(vTE)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 代入上边界u|AB=1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PTE=zeros(NPSUMTE,1);
for i=1:NXTE+1 %模型表层单元左上角点的编号
j=1+(i-1)*(NYTE+1);
vTE(j,j)=vTE(j,j)*10^10;
PTE(j)=vTE(j,j)*1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 稳定双共轭梯度法求解线性方程组 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tol=1e-15;%精度
maxit=100;%最大迭代次数
[LTE,UTE]=luinc(vTE,0);%不完全LU分解
xTE(:,ff)=bicgstab(vTE,PTE,tol,maxit,LTE,UTE);
end
size(xTE)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 计算视电阻率，张量阻抗及相位 %%%%%%%%%%%%%%%%%%%%%%%%%%
vect=(2+KONGQILAY+NYTE):NYTE+1:(NPSUMTE-NYTE+KONGQILAY-NYTE-1);%两列网格地表节点的编号
u1=zeros(size(vect,2),size(f,2));%改成列向量？
u2=zeros(size(vect,2),size(f,2));
u3=zeros(size(vect,2),size(f,2));
%u4=zeros(size(vect,2),size(f,2));
z=zeros(size(vect,2),size(f,2));
pc=zeros(size(vect,2),size(f,2));
phase=zeros(size(vect,2),size(f,2));
ve=[(1+KONGQILAY):NYTE:(NYTE*NXTE-NYTE+1+KONGQILAY),NYTE*NXTE-NYTE+1+KONGQILAY];%多加一个单元是为了补偿向量ve的长度与vect相同，方便后面计算
for j=1:size(f,2)
for i=1:size(vect,2)
u=(4e-7)*pi;
w=2*pi*f(j);
% u1(i,j)=xTE(vect(i),j);
% u2(i,j)=xTE(vect(i)+1,j);
% u3(i,j)=xTE(vect(i)+2,j);
% u4(i,j)=xTE(vect(i)+3,j);
u1(i,j)=xTE(vect(i),j);%
u2(i,j)=xTE(vect(i)-1,j);
u3(i,j)=xTE(vect(i)-2,j);
%u4(i,j)=xTE(vect(i)-3,j);
%四个节点的边长
%l=3*h;
l=2*a;
ux(i,j)=(u1(i,j)-u3(i,j))/l;
%(-11*u1(i,j)+18*u2(i,j)-9*u3(i,j)+2*u4(i,j))/(2*l); %近似求导数，对y求导
%z(i,j)=xTE(vect(i),j)*(sqrt(-1)*w*u)/ux(i,j);%阻抗
%pc(i,j)=abs(-sqrt(-1)*w*u*(xTE(vect(i),j)/ux(i,j))^2);%视电阻率
%phase(i,j)=atan(imag(z(i,j))/real(z(i,j)))*180/pi; %阻抗相位
HEY(i,j)=ux(i,j)/(sqrt(-1)*w*u);
%HY(i,j)=uy(i,j)/(sqrt(-1)*w*u);
end
end
plot(f,abs(HEY(30,:)));
%disD=0;%表水平距离
fid1=fopen('E:\workroom\无地形断面（TE，整体）.dat','wt+');
%disD=0;%表水平距离
for i=1:size(pc,1)
if i==1
disD(i)=0;
end
if i>1
disD(i)=disD(i-1)+a;
end
for j=1:size(f,2)
%fprintf(fid1,'\n%10.1f%10.1f%10.3f%10.2f%10.2f%10.2f',disD(i),i,log2(f(j)),pc(i,j), abs(phase(i,j)),abs(HEY(i,j)));%这个就是Hx
fprintf(fid1,'\n%10.1f%10.1f%10.3f%10.3f',disD(i),i,log2(f(j)),abs(HEY(i,j)));%这个就是Hx
end
end
% fid2=fopen('E:\workroom\无地形频点曲线（TE,整体）.dat','wt+');
% for j=1:size(f,2)
% %fprintf(fid2,'\n%10.0f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f',j,pc(j,1),pc(j,2),pc(j,3),pc(j,4),pc(j,5),pc(j,6),pc(j,7),pc(j,8),pc(j,9),pc(j,10),pc(j,11),pc(j,12),pc(j,13),pc(j,14),pc(j,15),pc(j,16),pc(j,17),pc(j,18),pc(j,19),pc(j,20),pc(j,1),pc(j,2),pc(j,3),pc(j,4),pc(j,5),pc(j,6),pc(j,7),pc(j,8),pc(j,9),pc(j,10),pc(j,11),pc(j,12),pc(j,13),pc(j,14),pc(j,15),pc(j,16),pc(j,17),pc(j,18),pc(j,19),pc(j,20));
% fprintf(fid2,'\n%10.0f%10.3f%10.3f%10.3f,%10.3f',j,f(j),pc(10,j),pc(20,j),pc(30,j));
% end
fclose('all');
评论0