geoid-height-correction.rar_geoid

clc %水准面不平行引起的高差改正 %已知高程(m) 点名：II405 高程：13.679 [names,lat] = textread('latitude.txt','%s %f'); %读入纬度概略值 %lat1 = (fix(lat)+fix((lat-fix(lat)*100)/60+(lat*100-fix(lat*100))/36)*pi/180; %将读入的dms纬度转化成弧度数 lat1 = fix(lat)*60+fix((lat-fix(lat))*100)+(lat*100-fix(lat*100))*100/60; [from,to,gaocha] = textread('height.txt','%s %s %f'); %读入高差观测值 单位 m n = size(names,1); % 观测点数目 m = size(from,1); % 观测高差数目 lev_no = zeros(m,2); % 每一观测值:起始点; 到达点 for k=1:m i = find(strcmp(names,from(k))~=0); % 起始点 j = find(strcmp(names,to(k))~=0); % 到达点 lev_no(k,1) = i; lev_no(k,2) = j; end A = zeros(m,n); M = zeros(m,n); H = zeros(n,1); H(2) = 13.679 ; % 高程初始值 eps = zeros(m,1); for k = 1:(m-1) i = lev_no(k,1); j = lev_no(k,2); if k ==1 H(i) = H(j)-gaocha(k); else H(j) = H(i)+gaocha(k); end end for k = 1:m i = lev_no(k,1); j = lev_no(k,2); if i~=0 A(k,i)=-1; end if j ~=0 A(k,j) = 1; end end for k = 1:m i = lev_no(k,1); j = lev_no(k,2); if i~=0 M(k,i)=1/2; end if j ~=0 M(k,j) = 1/2; end end delta_L = A*lat1; Hm = M*H; %两点平均高程 Lm = M*lat1; %两点平均纬度 for k = 1:n eps(k) = -0.0000015395*sin(2*((Lm(k)/60)*pi/180))*delta_L(k)*Hm(k)*1000; %单位 mm end disp('各点水准面不平行的高差改正eps为：单位 mm '); disp(eps); fid = fopen('result.txt','wt'); fprintf(fid,'水准面不平行的高差改正结果：'); fprintf(fid, '\n****************************************************'); fprintf(fid, '\n起点\t\t 终点\t\t改正值 单位：mm'); for k = 1:m fprintf(fid, '\n%s',cell2mat(from(k))); fprintf(fid,'\t\t'); fprintf(fid,'%s',cell2mat(to(k))); fprintf(fid,'\t\t'); fprintf(fid,'%6.4f',eps(k)); end fclose(fid);