学习记录——纵横断面面积计算
作者：互联网
最近在逐步学习MatLab，同时因为一些活动，要进行相关代码的编写，目前做了一个半成品出来，在这儿记录一下学习的过程吧。
代码中主要以大量的循环为主，因为对于MatLab的基础函数以及对于程序编写中参数类型的不熟悉，导致了很多的冗余运算，今后还要将努力学习。
clc
clear
close all
%% 读取所有点信息并展示
Indata=importdata('data.txt');
% fid=fopen('data.txt','r');
% data=fscanf(fid,'%s%d,%f,%f,%f\n');
Allpoint_data=Indata.data;
Azimuth_AB=Indata.textdata(3:4);
Point_name=Indata.textdata(5:end);
H0_Location=cell2mat(Indata.textdata(1));
A=isstrprop(H0_Location,'digit');
H0=str2num(H0_Location(A))/100000;
n=length(Allpoint_data); %存储数组总长，即点的总个数
X=Allpoint_data(:,1);
Y=Allpoint_data(:,2);
H=Allpoint_data(:,3);
plot(X,Y,'*y')
%% 提取K0/K1/K2的位置
SearchTool=strfind(Point_name,'K');
index=1;
for i=1:length(Allpoint_data)
if isempty(SearchTool{i}) == 0
key_point_Location(index,1)=i;
index=index+1;
end
end
%% 读取关键点的坐标
Key_point=zeros(length(key_point_Location),3);
for i=1:length(key_point_Location)
Key_point(i,:)=Allpoint_data(key_point_Location(i),:);
end
hold on
plot(Key_point(:,1),Key_point(:,2),'-b')
%% 计算纵断面
% 1.计算纵断面的长度
D=zeros(length(Key_point)-1,1);
for i=1:length(Key_point)-1
D(i)=sqrt((Key_point(i+1,1)-Key_point(i,1))^2+(Key_point(i+1,2)-Key_point(i,2))^2);
end
% 2.计算内插点的平面坐标
%（1）先计算内插点总个数
%（2）建立方位角矩阵Azimuth_Key_point，存储所有方位角信息
%（3）根据循环计算内插点的平面坐标
deltC=10;
%内插点个数
C_num=ceil(sum(D)/10);
%方位角矩阵
Azimuth_Key_point=zeros(length(Key_point)-1,1);
for i=1:length(Key_point)-1
%该处方位角计算是按横坐标为X，纵坐标为Y计算
deltX=Key_point(i+1,1)-Key_point(i,1);
deltY=Key_point(i+1,2)-Key_point(i,2);
if deltX>0 && deltY>0
Azimuth_Key_point(i)=atan(deltX/deltY);
else if (deltX>0 && deltY<0) || (deltX<0 && deltY<0)
Azimuth_Key_point(i)=atan(deltX/deltY)+pi;
else
Azimuth_Key_point(i)=atan(deltX/deltY)+2*pi;
end
end
end
Azimuth_Key_point=rad2deg(Azimuth_Key_point);
%计算内插点坐标
Index1=1;
for i=1:C_num
L=i*deltC;
%判别内插点所在位置
if L<D(Index1)
Interpolation_point(i,1)=Key_point(Index1,1)+L*sind(Azimuth_Key_point(Index1));
Interpolation_point(i,2)=Key_point(Index1,2)+L*cosd(Azimuth_Key_point(Index1));
else if L>D(Index1)
%每当进入新的Kj-Kj+1时，判别内插点是否在在关键点连线之外
if length(D)>=Index1+1
if L<sum(D(1:Index1+1))
D0=sqrt((Key_point(Index1+1,1)-Key_point(1,1))^2+(Key_point(Index1+1,2)-Key_point(1,2))^2);
Interpolation_point(i,1)=Key_point(Index1+1,1)+(L-D0)*sind(Azimuth_Key_point(Index1+1));
Interpolation_point(i,2)=Key_point(Index1+1,2)+(L-D0)*cosd(Azimuth_Key_point(Index1+1));
else
Index1=Index1+1;
end
else if length(D)<Index1+1
break
end
end
end
end
end
hold on
plot(Interpolation_point(:,1),Interpolation_point(:,2),'*k');
% 3.计算内插点高程
least_num=5;
for j=1:length(Interpolation_point)
for i=1:n
Dip(i,1)=sqrt((Interpolation_point(j,1)-X(i))^2+(Interpolation_point(j,2)-Y(i))^2);
end
[B,I]=sort(Dip);
I1=I(1:5);
Interpolation_point_height(j,1)=sum(H(I1(:))./B(1:5))/sum(1./B(1:5));
end
%4.计算纵断面面积
%（1）先将所有关键点放入内插点矩阵中的对应位置
%（2）按面积公式计算
Dc=D./deltC;
Dc=floor(Dc);
Vertical_section_point=zeros(length(Key_point)+length(Interpolation_point),3);
Vertical_section_point(1,1)=Key_point(1,1);
Vertical_section_point(1,2)=Key_point(1,2);
Vertical_section_point(1,3)=Key_point(1,3);
Index2=1;
Index3=1;
for i=2:length(Vertical_section_point)
if i~=sum(Dc(1:Index2))+1+Index2
Vertical_section_point(i,1)=Interpolation_point(Index3,1);
Vertical_section_point(i,2)=Interpolation_point(Index3,2);
Vertical_section_point(i,3)=Interpolation_point_height(Index3,1);
Index3=Index3+1;
else
if Index2<=length(Dc)
Index2=Index2+1;
Vertical_section_point(i,1)=Key_point(Index2,1);
Vertical_section_point(i,2)=Key_point(Index2,2);
Vertical_section_point(i,3)=Key_point(Index2,3);
else
break
end
end
end
% 每个梯形计算面积
S=zeros(length(Vertical_section_point)-1,1);
for i=1:length(Vertical_section_point)-1
deltL=sqrt((Vertical_section_point(i,1)-Vertical_section_point(i+1,1))^2+(Vertical_section_point(i,2)-Vertical_section_point(i+1,2))^2);
S(i,1)=deltL*(Vertical_section_point(i,3)+Vertical_section_point(i+1,3)-2*H0)/2;
end
%计算总面积
OutS=sum(S);
%% 计算横断面
% 1.计算横断面中心点
Central_point=zeros(length(Key_point)-1,2);
for i=1:length(Key_point)-1
Central_point(i,1)=(Key_point(i,1)+Key_point(i+1,1))/2;
Central_point(i,2)=(Key_point(i,2)+Key_point(i+1,2))/2;
end
hold on
plot(Central_point(:,1),Central_point(:,2),'*r');
%求得每根过中点垂直于关键点连线的直线斜率
Mid_K=zeros(length(Key_point)-1,1);
B=zeros(length(Key_point)-1,1);
for i=1:length(Key_point)-1
Mid_K(i)=-1*(Key_point(i+1,1)-Key_point(i,1))/(Key_point(i+1,2)-Key_point(i,2));
B(i)=Central_point(i,2)-Central_point(i,1)*Mid_K(i);
end
%计算内插点坐标
deltCS=5;
Interpolation_point_CS=zeros(ceil(25*2/deltCS*length(Central_point)),3);
%(1)根据斜率计算方位角
Azimuth_Central_point=zeros(length(Mid_K),1);
for i=1:length(Central_point)
Azimuth_Central_point(i)=atan(Mid_K(i));
end
% Azimuth_Central_point=rad2deg( Azimuth_Central_point);
%（2）计算两边的内插点坐标
for j=1:length(Azimuth_Central_point)
K_par=(2*j-2)*25/deltCS;
K_par1=(2*j-1)*25/deltCS;
Index4=1;
for i=1:deltCS
%判断方位角所在象限
if Azimuth_Central_point(j)>0 && Azimuth_Central_point(j)<pi/2
Interpolation_point_CS(K_par+Index4,1)=Central_point(j,1)+i*deltCS*cos(Azimuth_Central_point(j));
Interpolation_point_CS(K_par+Index4,2)=Central_point(j,2)+i*deltCS*sin(Azimuth_Central_point(j));
Interpolation_point_CS(K_par1+Index4,1)=Central_point(j,1)-i*deltCS*cos(Azimuth_Central_point(j));
Interpolation_point_CS(K_par1+Index4,2)=Central_point(j,2)-i*deltCS*sin(Azimuth_Central_point(j));
else if Azimuth_Central_point(j)>pi/2 && Azimuth_Central_point(j)<pi
Interpolation_point_CS(K_par+Index4,1)=Central_point(j,1)+i*deltCS*cos(Azimuth_Central_point(j));
Interpolation_point_CS(K_par+Index4,2)=Central_point(j,2)-i*deltCS*sin(Azimuth_Central_point(j));
Interpolation_point_CS(K_par1+Index4,1)=Central_point(j,1)-i*deltCS*cos(Azimuth_Central_point(j));
Interpolation_point_CS(K_par1+Index4,2)=Central_point(j,2)+i*deltCS*sin(Azimuth_Central_point(j));
else if Azimuth_Central_point(j)>pi && Azimuth_Central_point(j)<3*pi/2
Interpolation_point_CS(K_par+Index4,1)=Central_point(j,1)-i*deltCS*cos(Azimuth_Central_point(j));
Interpolation_point_CS(K_par+Index4,2)=Central_point(j,2)-i*deltCS*sin(Azimuth_Central_point(j));
Interpolation_point_CS(K_par1+Index4,1)=Central_point(j,1)+i*deltCS*cos(Azimuth_Central_point(j));
Interpolation_point_CS(K_par1+Index4,2)=Central_point(j,2)+i*deltCS*sin(Azimuth_Central_point(j));
else if Azimuth_Central_point(j)>pi/2 && Azimuth_Central_point(j)<pi
Interpolation_point_CS(K_par+Index4,1)=Central_point(j,1)-i*deltCS*cos(Azimuth_Central_point(j));
Interpolation_point_CS(K_par+Index4,2)=Central_point(j,2)+i*deltCS*sin(Azimuth_Central_point(j));
Interpolation_point_CS(K_par1+Index4,1)=Central_point(j,1)+i*deltCS*cos(Azimuth_Central_point(j));
Interpolation_point_CS(K_par1+Index4,2)=Central_point(j,2)-i*deltCS*sin(Azimuth_Central_point(j));
end
end
end
end
Index4=Index4+1;
end
end
hold on
plot(Interpolation_point_CS(1:10,1),Interpolation_point_CS(1:10,2),'-b')
plot(Interpolation_point_CS(11:20,1),Interpolation_point_CS(11:20,2),'-b')
% plot(Interpolation_point_CS([1,6],1),Interpolation_point_CS([1,6],2),'*k')
% 将所有横断面上的点依次放入矩阵Cross_section_point
Cross_section_point=zeros(length(Interpolation_point_CS)+length(Central_point),3);
Index5=1;
Index6=1;
for i=1:length(Cross_section_point)
K_par2=25/deltCS+1+(Index6-1)*(25/deltCS*2+1);
if i~= K_par2
Cross_section_point(i,1)=Interpolation_point_CS(Index5,1);
Cross_section_point(i,2)=Interpolation_point_CS(Index5,2);
Cross_section_point(i,3)=Interpolation_point_CS(Index5,3);
Index5=Index5+1;
else
Cross_section_point(i,1)=Central_point(Index6,1);
Cross_section_point(i,2)=Central_point(Index6,2);
Cross_section_point(i,3)=0;
Index6=Index6+1;
end
end
hold on
plot(Cross_section_point(:,1),Cross_section_point(:,2),'*k')
%计算内插点高程
for j=1:length(Cross_section_point)
for i=1:n
Dip_CS(i,1)=sqrt((Cross_section_point(j,1)-X(i))^2+(Cross_section_point(j,2)-Y(i))^2);
end
[B_CS,I_CS]=sort(Dip_CS);
I1_CS=I_CS(1:5);
Cross_section_point(j,3)=sum(H(I1_CS(:))./B_CS(1:5))/sum(1./B_CS(1:5));
end
%计算横断面面积
S_CS=zeros((length(Cross_section_point)-length(Central_point))/2,length(Central_point));
Index8=1;
for i=1:length(Central_point)
Index7=1;
for j=1:(length(Cross_section_point)-length(Central_point))/2
K_par3=(i-1)*(25/deltCS*2+1);
deltD=sqrt((Cross_section_point(j+ K_par3,1)-Cross_section_point(j+ K_par3+1,1))^2+(Cross_section_point(j+ K_par3,2)-Cross_section_point(j+ K_par3+1,2))^2);
S_CS(Index7,Index8)=(Cross_section_point(j+ K_par3,3)+Cross_section_point(j+ K_par3+1,3)-2*H0)/2*deltD;
Index7=Index7+1;
end
Index8=Index8+1;
end
S_M0=sum(S_CS(:,1));
S_M1=sum(S_CS(:,2));
标签：point,记录,学习,length,Location,Key,Azimuth,横断面,data 来源： https://blog.csdn.net/m0_49646774/article/details/119011306