赞
踩
% smt2djxfwd01.m
% 二维大地电磁有限元数值模拟 矩形+线性插值
clear
close all
clc
format short g
%
% 读频率数据
%%
[nfreq,freq]=sinfreq('freq.dat');
%
% 读模型数据
%
[idx,ney,nez,dy,dz,res,resbottom]=sinmt2dmodjx1('mod2dfem1')
t1=cputime;
for ifreq=1:nfreq
ifreq;
frek=freq(ifreq);
[gk,s]=sglobalksrectlinear(frek,idx,ney,nez,dy,dz,res,resbottom);
u=gk\s;
%
% 计算视电阻率
%
[es,hs,ksair] = sfldsurf(idx,frek,u,ney,nez,dy,dz,res);
[z,r,p]=sImpsurf(idx,frek,es,hs,ney);
zimp(ifreq,1:ney+1) = z(1:ney+1,1);
rhoa(ifreq,1:ney+1) = r(1:ney+1,1);
phaz(ifreq,1:ney+1) = p(1:ney+1,1);
end
t2=cputime-t1
stemp = sprintf('花费时间: %g s',t2);
disp(stemp);
clear x
x(1)=0;
for iy=2:ney+1
x(iy) = x(iy-1)+dy(iy-1);
end
for iy=1:ney+1
x(iy) = x(iy) - x(ney+1)/2;
end
clear x1
for i=6:16
x1(i-6+1) = x(i)/1000;
end
clear x
x = x1;
nx = length(x);
clear y
for i=1:nfreq
y(i) = log10(freq(i));
end
clear x1 y1 z1
for i=1:nx
for j=1:nfreq
x1(i,j) = x(i);
y1(i,j) = y(j);
z1(i,j) = log10(rhoa(j,i-1+6));
end
end
figure(1);
colormap(jet);
contourf(x1,y1,z1);
title('TM模式的视电阻率拟断面图');
xlabel( ' x /km');
ylabel(' log10(f) /Hz');
colorbar('vertical');
clear y
y(1)=0;
for i=2:nez+1
y(i) = y(i-1) + dz(i-1);
end
for i=1:nez+1
y(i) = y(i)/1000;
end
clear x1 y1 z1;
for i=1:nx
for j=1:nez-6
x1(i,j) = x(i);
y1(i,j) = -y(j);
z1(i,j) = res(j,i-1+6);
end
end
figure(2);
colormap(jet);
pcolor(x1,y1,z1);
title('TM模式网格剖分示意图');
xlabel(' x/km');
ylabel(' z/km');
axis equal;
axis tight;
% colorbar('vertical');
D140
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。