WDEM接口小程序

fid1 = fopen('xy.dat','rt');   
fid2 = fopen('elevent0.dat','rt');
fgetl(fid1);
b = fscanf(fid1,'%f %f %f %f %f %f %f %f %f %f %f %f %f %f',[14,inf]);
Ele = fscanf(fid2,'%f %f',[2,inf]);
fclose(fid1);fclose(fid2);
b = b';
Ele = Ele';
%% 文件头 根据实际情况改动
Num_Edge = 16;     %扩边点数
Coe_Edge = 1.3;    %扩边系数
Terrain = 1;       %0：无地表地形数据，1：有地表地形数据
Err = 0           ;%0：无误差数据，1：有误差数据
source = 1;        %0：数据无源，1：数据有源，
site_source = [3269728.001 756030.001 1264.663    %场源A（x,y,z)
               3270773.001 757175.001 1266.558];  %场源B（x,y,z)
%% 判断频点数
if abs( sum(b(:,1))/length(b)-b(1,1)) < 1e-4
    n_fre=length(b);
else
    ik=2;
    while abs(b(ik,1)-b(1,1) )< 1e-2;
        ik=ik+1;
    end
    n_fre=ik-1;
end
%%                        
n_site = length(b)/n_fre;         %测点数
r_MN = abs(b(1+n_fre,1)-b(1,1));  %点距
site_x = (b(:,6)+b(:,8))/2;       %测点X坐标
site_y = (b(:,7)+b(:,9))/2;       %测点y坐标
j = 1;
m = 1;
site_filexy = b(:,1);

for i=1:n_fre:length(b)
    site_x0(j) = site_x(i);
    site_y0(j) = site_y(i);
    site(j) = site_filexy(i);
    j=j+1;
end
if length(site)==length(Ele)
else
    for k=1:length(site)
        gaocheng(m)= Ele((find(Ele(:,1)==site(k))),2);
        m = m+1;
    end
end

site_xyz = [site_x0',site_y0',gaocheng'];%
MainData = [b(:,1),b(:,10),b(:,11),b(:,12)];%/(r_MN/1000)
FileTitle = [n_site,r_MN,Num_Edge,Coe_Edge,n_fre,Terrain,Err]';


filename = 'WDEM_input.dat';
dlmwrite(filename,FileTitle, 'delimiter', ' ');
dlmwrite(filename,site_xyz, '-append', 'roffset', 3, 'delimiter', ' ', 'precision', '%f');
dlmwrite(filename,MainData, '-append', 'roffset', 1, 'delimiter', ' ', 'precision', '%.6f');
dlmwrite(filename,source, '-append', 'roffset', 1, 'delimiter', ' ');
if source == 1
dlmwrite(filename,site_source, '-append', 'precision', '%.4f', 'delimiter', ' ');  
end