tdu70.m

重力异常的向上延拓程序

%二度异常向上延拓(利用剖面上的异常值)matlab7.4计算程序
%Design by 赵玉岩
%最后修改日期：2008年4月21日23:03:19
%Δg原始数据为excel文件，内容格式为：
%点号         Δg   
%-2           0.00
%-1           0.00
%0            0.00
%1            0.00
%2            0.00
%Bi原始数据为，内excel文件容格式为：
%Bi
%0.00
%0.00
%0.00
%##########################################################################
clear;
clc;%清空窗体和内存内容
%以下读入相关参数的对话框
%以下是读入点个数的对话框
cydj=input('输入采样点距')                      %采样点距
ytgd=input('输入延拓高度')                      %延拓高度
if ytgd<cydj 
    msgbox '延拓高度小于采样点距，程序无法工作'   
    break
end    
if rem(ytgd,cydj)~=0
    msgbox '延拓高度不是采样点距的整数倍，程序无法工作'   
    break
end 
[datafile,pathname]=uigetfile({'*.xls';'*.*'},'打开Δg原始数据文件：');%打开Δg原始数据文件
[ydata headtext]= xlsread([pathname datafile]);
subplot(2,1,1),plot(ydata(:,1),ydata(:,2))
hold on
ydx=size(ydata);%矩阵大小
yhs=ydx(1);%矩阵行数
ydh=ydata(:,1);
ydg=ydata(:,2);
[datafile,pathname]=uigetfile({'*.xls';'*.*'},'打开系数原始数据文件：');%打开系数原始数据文件
[bdata headtext]= xlsread([pathname datafile]);
bdx=size(bdata);%矩阵大小
bhs=bdx(1)-1;%矩阵行数
if 2*bhs>yhs
    msgbox '系数值个数大于测点的1/2，程序无法工作'   
    break
end   
%对系数进行处理
bn=[];
for i=0:bhs-1
    bn=[bn bdata(bhs-i)];
end
bn=[bn bdata'];    
%ydg进行扩边处理
nydg=td(ydg',ytgd/cydj*bhs);
ydg=nydg';
yhs=yhs+ytgd/cydj*bhs*2;
%以下开始计算,不包括边缘数据
ytdgjz=[];
for i=(ytgd/cydj*bhs+1):(yhs-ytgd/cydj*bhs)
    ndg=[];
   for j=(i-ytgd/cydj*bhs):(ytgd/cydj):(i+ytgd/cydj*bhs)
       ndg=[ndg ydg(j)];
   end
   ytdg=ndg*(bn');
   ytdgjz=[ytdgjz ytdg];
end
[filename2, pathname2] = uiputfile( ...
       {'*.xls', '电子表格(*.xls)'}, ...
        '二度异常向上延拓计算结果：');
a1={'点号' '向上延拓值'};  
xlswrite([pathname2 filename2],a1,'sheet1','A1');   
xlswrite([pathname2 filename2], ydh ,'sheet1','A2');   
xlswrite([pathname2 filename2], ytdgjz' ,'sheet1','B2');
subplot(2,1,2),plot(ytdgjz)
msgbox '计算完成！'
hold off