赤道子午线弧长反演大地纬度

最新推荐文章于 2024-09-28 22:35:00 发布

黑猫巫师

最新推荐文章于 2024-09-28 22:35:00 发布

阅读量753

点赞数 5

分类专栏：大地测量基础文章标签： matlab

本文链接：https://blog.csdn.net/weixin_45635194/article/details/109134249

版权

大地测量基础专栏收录该内容

1 篇文章 0 订阅

订阅专栏

matlab赤道子午线弧长反演大地纬度

话不多说，直接看公式
太懒了，直接拍课本的照片
在这里插入图片描述

在这里插入图片描述

clear
clc
S=input('请输入由赤道起算的子午线弧长S:');
a=input('请输入长半轴半径a：');
f=input('请输入扁率f：');
e2=2*f-f*f;
[A,B,C,D,E,F,G]=T(e2);     %调用T函数产生运算系数
n=0;
b2=S/(a*(1-e2)*A);
s12=a*(1-e2)*(A*(b2)-B*sin(b2)*cos(b2)-C*(sin(b2))^3*cos(b2)-D*(sin(b2))^5*cos(b2)-E*(sin(b2))^7*cos(b2)-F*(sin(b2))^9*cos(b2)-G*(sin(b2))^11*cos(b2));
while abs(S-s12)>0.001
    n=n+1;
    b2=b2+(S-s12)/(a*(1-e2)*(A-B*cos(2*b2))-C*(sin(b2))^2*(cos(2*b2)+2*(cos(b2))^2)-D*(sin(b2))^4*(cos(2*b2)+4*(cos(b2))^2)-E*(sin(b2))^6*(cos(2*b2)+6*(cos(b2))^2)-F*(sin(b2))^8*(cos(2*b2)+8*(cos(b2))^2));
    s12=a*(1-e2)*(A*(b2)-B*sin(b2)*cos(b2)-C*(sin(b2))^3*cos(b2)-D*(sin(b2))^5*cos(b2)-E*(sin(b2))^7*cos(b2)-F*(sin(b2))^9*cos(b2)-G*(sin(b2))^11*cos(b2));
  end    
B2=r2d(b2);
fprintf('\n所求大地纬度B=%.9f\n\n',B2);


function [A,B,C,D,E,F,G]=T(e2)
A=1+3/4*e2+45/64*e2^2+175/256*e2^3+11025/16384*e2^4+43659/65536*e2^5+693693/1048576*e2^6;
B=3/4*e2+45/64*e2^2+175/256*e2^3+11025/16384*e2^4+43659/65536*e2^5693693/1048576*e2^6;
C=15/32*e2^2+175/384*e2^3+3675/8192*e2^4+14553/32768*e2^5+231231/524288*e2^6;
D=35/96*e2^3+735/2048*e2^4+14553/40960*e2^5+231231/655360*e2^6;
E=315/1024*e2^4+6237/20480*e2^5+99099/327680*e2^6;
F=693/2560*e2^5+11011/40960*e2^6;
G=1001/4096*e2^6;
end

   function deg=r2d(rad)
    %   弧度转角度
    rad1=rad*180/pi;
    deg1=fix(rad1);
    min=fix((rad1-deg1)*60);
    sec=((rad1-deg1)*60-min)*60;
    deg=deg1+min/100+sec/10000;
   end